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sheath. The non-existence of a one-dimensional sheath 
solution is proved and a computer model with 
two-dimensional, periodic active sites representing a flat 
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the effects of a magnetic field and Joule heating are 
studied. Results are compared with experimental 
observations. To supplement the sheath investigation a 
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I. 



INTRODUCTION 



A. BACKGROUND 

A Magnetchydrody namic (MHD) generator is a device for 
converting the energy of a flowing gas, or of a liquid 
metal, directly into electrical energy. As with 
conventional electrical generators, the conducting medium 
crosses a magnetic field; however, MHD deals with conducting 
fluids instead of solid moving parts. Anticipated 
efficiencies are far above those of conventional generators 
because of the higher temperatures in the conversion 
channel. The thermal energy of flowing gases is changed 
directly to electrical energy, eliminating the mechanical 
energy step cf the conventional generator which means that 
there are fewer parts to wear out. 

Although MHD power generation has been studied for 
nearly 30 years, financing of HUD research waned in the 
1960's. Interest is now being rekindled because of a new 
"energy consciousness" of government, industry, and the 
military. MHD offers an attractive alternative energy 
transformation means which is relevant to the pursuit of 
clean ways of using the vast coal reserves in the United 
States and aaking more efficient use of dwindling petroleum 
resources. MHD promises to fulfill the requirements for a 
light-weight, high-power source for military use. For 
example, the Navy, recognizing its energy requirements and 
dependency on world fossil-fuel supplies[ 1 ], is sponsoring 
research on practical MHD devices with large power outputs 
and high eff iciencies[ 2 ]. The Air Force is interested in 
power sources for airborn weapons systems. 

The basic MHD device consists of a channel through which 
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ionized gas or liquid metal flows. A magnetic field is 
imposed perpendicular to the flow direction; as a result, an 
electric field is induced and pov/er is tapped by means of 
electrodes. Principal concern is with ionized gas 
(hereafter refered to as plasma) devices since with these 
the combustion energy can be more directly changed to 
electrical power. 

A major problem involved in the performance of MHD 
generators is the high voltage loss in the vicinity of the 
electrodes. Accurate performance predictions in the design 
of HHD devices require a realistic determination of these 
voltage losses. This research investigates the nature of 
electrode voltage drops and computer models are developed to 
describe the physical phenomena involved. With the results 
from this investigation, appropriate steps may be taken to 
minimize the losses. 

The principal voltage loss mechanisms in the KHD 
generator can be divided into two main classes, ohmic and 
sheath losses. Ohmic drops are those that occur because of 
the finite conductivity of a real plasma. Thermal boundary 
layers, degree and kinetics of ionization, and Joule heating 
are factors affecting the ohmic resistivity of the plasma. 
Sheath drops occur as the result of the Debye shielding 
which forms a non-neutral layer adjacent to the electrode 
and results in a space charge field. ilaterial problems 
restrict the temperatures at which the electrodes can 
operate. In many cases cooling of the electrodes is 
required since the plasma, in order to maintain a high 
ionization, must be hotter than the working temperatures of 
most materials. This temperature difference between the 
electrodes and the plasma further aggravates the voltage 
losses because of the presence of the thermal boundary 
layers. As will be shown, voltages losses can be as much as 
50% or more of the total power output, with sheath drops 
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accounting for a non-negligible fraction of the drop and 
boundary layer losses the rest. 

B. REVIEW Of PREVIOUS WORK 

Analyses and experiments on boundary layer phenomena in 
MHD generators are numerous. It is not the purpose of this 
work to further add to the literature on the boundary layer 
losses, but, as will be explained in Part C of this section, 
this phase of the problem was treated as a necessary 
appendage to the problem of the sheath drop. Consequently a 
short review cf boundary layer work is presented here: High 

and Felderman[3] as well as Argy ropoulos, Deraetriades, and 
Lackner[4] studied the problem of the turbulent boundary 
layer and several descriptions of the nature of the 
temperature field. Doss, Dwyer and Hoffman[5] investigated 
the boundary layer and used a rudimentary collisionless 
sheath as an '‘inside” boundary condition. Kessler and 
Eustis[6] reported on experiments with electrode temperature 
effects on turbulent boundary layers, and Rubin and 
Eustis[7] extended Kessler's work to include the effects of 
electrode size. Wu et al.[3] and Oliver and Hitchner[9] 
investigated non- uniformities of current distributions 
within the boundary layer. 

Although the existence of the sheath is well understood, 
its effects have been investigated to a much lesser extent 
than those of the boundary layer. This is principally 
because this loss is described by a relatively complicated 
set of coupled, non-linear, partial differential equations 
that present considerable difficulty for' numerical 
solutions. There is no known exact solution for these 
equations. 

Most work on sheath phenomena is embodied in "probe” 
theory i! . estigations, that is in the mutual effects on a 
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quiescent plasma which is disturbed by the presence of an 
electric prcbe. Such work is relevant to MHD electrodes 

since the anode is essentially a heavily biased probe in 
contact with a plasma, which is quiescent within the sheath 
reqion (see Part C of this section) . The cathode, which is 
an electron emitter, has some of the same characteristics as 
the anode, but is more complicated to represent 
analytically. The description of the cathode will not be 
attempted here. 

Necessary simplifications limit existing probe solutions 

to special cases as can be seen by sampling the literature 

regarding plasma probes. Lam[10] solves a sheath problem by 

matching the boundaries of "inner" and "outer" solutions and 

discusses their dimensionality. Stahl and 3u[11] use the 

same approach of separating the sheath, ambipolar region and 

free stream. Additionally, they prove the existence of a 

sheath on a flat probe. Cohen[12] criticizes this approach 

of separate regions on mathematical grounds saying that 

quantities are forced to fit and some derivatives tend to be 

discontinuous. McKee and Hitchner[13] deal with a 

collisionless sheath ) but include ionization and 

s 

recombination in the ambipolar region. Bailey and 
Touryan[14] investigate a sheath that is large enough to be 
of the same order of magnitude as the boundary layer, and 
take advantage of Blasius' similitude co-ordinates to reduce 
the problem to one dimension. 

Several solutions are available for spherical probes 
including Kiel[15J , Barad and Cohen[16] , Su and Lam[17] , 
and Cohen[12]. Kiel ignored diffusion while Barad and Cohen 
neglected ionization and recombination. Su and Lam, and 
Cohen were the first to use a systematic analysis of probes 
in collision-dominated plasmas. Lengyel[18], by dropping 
diffusion, was able to modify the elliptic equations in such 
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a way as to solve them by the method of characteristics as 
one might do with hyperbolic equations. 

Chung/ Talbot, and Touryan[19], in a rather extensive 
review of probe work, state that no general solution is 
available for determining charge density and species 
temperature for probes small relative to boundary layer 
thicknesses. The work of this thesis may help to fill that 
void because the model of the anode is that of an array of 
point or line probes immersed in the non-convective portion 
of a boundary layer. Results of this work nay also be 
useful to other fields such as arc-discharges and lasers. 



C. PHYSICAL MODELING 

Both theoretical and experimental results indicate that 
in the HHD environment, where pressures are near 
atmospheric, the sheath lies well within the boundary layer 
(about 10 ~s Hi for the sheath thickness or three orders of 
magnitude less than the boundary layer) . Additionally, the 
sheath is only about one or two orders of magnitude larger 
than the electron mean-f ree-path for combustion gases. It 
becomes apparent then, that convection plays little or no 
role in the behavior of the sheath other than to modify the 
gross temperatures encountered near the wall. Consequently, 
the boundary layer or ohmic problem can be divorced from the 
sheath and treated separately[ 20 ]. The boundary layer, as 
well as the wall conditions, determines the gas temperature 
for the sheath region. Thus, the boundary layer voltage 
drop can be added to the sheath drop to give the total loss 
due to electrode effects. 

As mentioned earlier, analyses of boundary layer effects 
in HHD have teen developed, but they tend to be complicated 
and difficult to use. In order to have a complete picture 
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of electrode drops and to assist in interpreting data, a 
simplified method for determining the voltage drop across a 
thermal boundary layer has been devised. It has the 
capability cf rendering quick results with little computer 
space and time. It does not pretend to be a complete 
analysis, but through appropriate assumptions and 
simplifications it gives the drop as a function of 
temperature only. Appendix A gives the details of the 
method along with its limitations. 

An investigation of the nature of the controlling 
equations for the sheath shows that a one-dimensional 
solution does not exist except for very special cases. 
This investigation resulted in the publication of Sef. 21, 
and a summary of that paper may be found in Appendix E. It 
provides a basis for many of the assumptions used in this 
work. The following conclusions which appear in Eef. 21 
helped to shape the model chosen for this work: Current 
density must decrease away from the electrode by whatever 
means possible. The mechanisms of geometry, current 
constrictions at the electrode, or ionizaticn/recombination 
can effect this decrease singly or jointly. A 
two-dimensional flat plate in the absence of 
ionization/recombination would be at least partially 
one-dimensional away from the edges and is therefore not a 
workable model. 

Since previous solutions have been obtained with 
spherical probes, it appears that a geometrical decrease of 
current density is a proper means of satisfying the nature 
of the equations. Chen[22] used ionization/recombination to 
extract a strictly one-dimensional solution for the sheath. 
What remains then is to study the features that result from 
current constrictions at discrete periodic sites rather than 
uniforia current density across the electrode surface. 
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D. 



OBJECTIVES OF THE PRESENT WORK 



In this work, a two-dimensional, flat plate, 
non-emitting electrode is modeled. Here, the current has 
only discrete points through which to flow on the electrode 
surface (See Fig. 1) . These points are assumed to be 
periodic. The points simulate the current constriction at 
the electrode which may be there because of temperature 
irregularities or surface roughness. A collection of these 
points or nodes might represent the anode spots observed by 
Biblarz[23] and Kinblin[24]. In a two-di mensional 
representation the point is actually a '’wire'* of unit depth. 
The spacing of the active sites must be greater than the 
sheath thickness in order to avoid one-dimensional effects, 
and small enough to be compatible with a two-dimensional 
computational field. It must be emphasized here that this 
work does not consider the phenomenon of arc discharge, 
which is the result of thermal instabilites[ 25 ]. Rather it 
looks at the current constriction required to satisfy 
continuity of charge. Ohm’s law, and Poisson's equation. 

The primary objective of this work is to present a 
mathematical model of a two-dimensional (Cartesian geometry) 
electrode in a quiescent HHD plasma. The controlling 
equations are solved for the potential, the current, and the 
charge density distributions in tne vicinity of the 
electrode. The additional effects of a magnetic field and 
Joule heating are investigated. It appears that this is the 
first investigation of the magnetic field within the sheath. 
The resulting program generates sheath and ambipolar regions 
in a self-consistent way using the same set of equations 
throughout the field, obviating the need to match layers. 
The size of the sheath and the voltage drop attributable to 
its shielding effects are investigated. 
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E. ORGANIZATION 



In Section II the physical concepts and the controlling 
and boundary equations are developed. The equations are 
cast into a form for use in the computer in Section III. It 
includes a description of the numerical methods v;hich 
produced successful results, as well as a short chronology 
of techniques that failed to work. A summary of t.he results 
of potential and charge distributions, and current and 
temperature plots for all the cases studied is included in 
Section IV. Section V lists the conclusions reached from 
both a physical and numerical point of view and gives 
examples on the use of the results. 

Appendix A shows the technique for determining the 
boundary layer losses. Appendix 3 is a condensation of 
Ref. 21 proving the non-existence of a one-dimensional 
solution for the limit of the non-reacting flux of charges 
in a collision-dominated plasma. A dimensional analysis 
which derives a set of independent non-dimensional variables 
for the controlling equations is found in Appendix C. 
Additionally, Appendix C contains a fractional analysis of 
these and other equations to estimate the significance of 
certain physical effects. Finally, the computer programs 
used in this work are included in the section "Computer 
Programs", with explanations as to their use. 
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II. ANA LYS IS OF THE SHEATH AND AHB IPOLAR REGIONS 
A. THE SHEATH ENVIRONMENT 

A sheath is a non-neutral region which lies adjacent to 
the electrode or insulator surface. The strict definition 
of a plasma requires charge neutrality, thus the sheath is 
not part of the plasma since a space charge exists. For 
example, the anode is positive, which means that it attracts 
electrons and repels positive ions. As the electrons 
collect at the anode surface, the anode potential is 
partially shielded from the rest of the plasma by the 
space-charge potential drop. 

Betv;een the sheath and the free stream lies the 
ambipolar region. As the name implies, it has an equal 
number of positive and negative charges. In this transition 
region the electrons are slowed by heavier ions, and strong 
concentration gradients exist. The electrodes themselves 
are usually metallic, and may be pins, wires, or flat plates 
separated by ceramic insulator wall segments (segmented 
electrodes) . The surface of the electrodes, even when highly 
polished have roughnesses of the order of the sheath 
thickness, and in addition exchange material with the 
flowing plasma, further adding to the irregularities. These 
irregularities increase the tendency for the current to 
constrict by providing active sites for the current to flow 
along minimum energy paths. 

The size of the sheath in tlHD generators is about 
10 -s ra, putting it well within the boundary layer whicn is 
about 10“2 m. The flow is therefore essentially stagnated 
in the sheath region. Though the ambipolar region extends 
further into the boundary layer, a fractional analysis 
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(Appendix C) on convection shows that the characteristic 
length for convection in most flows is much greater than the 
conduction or diffusion lengths within the sheath and 
ambipolar regions. Consequently, convection may be entirely 
assigned to the ohmic study of electrode losses (see 
Appendix A) . The boundary layer lies within the undisturbed 
plasma. That is, it lies beyond the ambipolar region, there 
is charge neutrality, and electron/ion concentrations are as 
would be predicted by the Saha equation[26] for plasmas at 
local equilibrium. 

In the model presented here, the following assumptions 
are made: 

1) Steady state, 

2) Chemical equilibrium in the boundary layer, but frozen 

flow in the sheath, 

3) No induced magnetic field, i.e., low Magnetic Reynolds 

number, 

4) Negligible ion slip, 

5) J << J , 

i e 

6) No continuum radiation losses in the energy equation, 

7) No ion emission, and 

8) Neutral and ion particle temperature is T , unaffected 

o 

by Joule heating and uniform within the domain of the 

sheath . 



B. CONTROLLING EQUATIONS 
1 . B^ic Equa tio ns 

For some MHD plasmas the mean-f ree-path is small 
enough that the sheath can be treated as a continuum. A set 
of collisional equations is then used to describe the 
sheath, ambipolar, and adjacent free stream regions. No 
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matching between these regions is necessary since they are 
generated in a self-consistent way within the computational 
array. The equations used are a combination of the 
conservation equations and Maxwell's equations. 
Specifically, they are expressions of conservation of 
charge, momentum and energy, as well as Gauss' law. 



A special representation of the electron momentum 
equation is embodied in a generalized Ohm's law and may be 
written [Ref. 26, p. 363] 

3 . 



J = -u n eVd) - 
o s ^ 



w 



(J„xB) ± D eVn 



. ( 1 ) 



where the + sign applies to the electrons and - is tor ions. 

The first terra on the right-hand-side represents conduction, 

the second is due to the magnetic JXB current and the third 

terra represents diffusion. The terra (3 is the Hall 

s 

parameter and gives the relative effect of the magnetic 
field in the plasma. The equation for species continuity 
can be written as 



V • J = n e 
s s 



( 2 ) 



Appendix C shows that 


the characteristic len 


gt h 


ioniz at ion/reconbi nation 


is much 


larger 


than the 


sheath 


ambipolar region lengths. 


Under 


these 


condition 


S Eq. 


becomes 











for 



V- J = 0 
s 



{ 3 ) 



Poisson's equation comes from Gauss* law for the divergence 
of an electric field and is given in potential form as 



V <p = -(n.-n^) e/c^ 



('0 



Another equation that will be useful later and which 
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applies under local equilibrium is the Einstein relation. 



D /y = kT /e 
s s s 



( 5 ) 



Electrons absorb energy at the rate of J *E per 

e 

unit volume from the electric field. Because of their 
relatively small mass, the electrons are inefficient in 
transfering energy to the other particles. For that reason 
much of this energy goes into raising the electron 
temperature above that of the heavier particles. Collisions 
with these particles, both elastic and inelastic, tend to 
limit the energy retained by the electrons. In a quiescent 
plasma an expression accounting for the electron energy 
balance is (Eef. 26, p. 240-243) 

J -E = 3n k(T -T) y 
e e e ^ 

s 



V. 



es 



6 m 
s e 

m 



( 6 ) 



2 • Hon -D im ens io n al Parameter s 

The non-dimensional parameters developed in Appendi: 
C are repeated here for use in this section. 



"i = 



V = LV 



n = n /n 
e e o 



X = x/L 



<P = <P/<P, 



y = y/L 



(7) 



v;here : 



(f)^ = kT^/e L = e^/(e^kT^) and = (kT^e^/e^) ^ 



T is the temperature of the neutral gas which is constant 
o 

throughout the region of interest, and independent of Joule 
heating. 
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Each of the three characteristic parameters , L 

o 

and n are functions only of a characteristic energy kT . 
o o 

ft now represents a non-dimensional electron charge, or a 
e 

measure of the deviation from the equilibrium density value 
in the free stream. ft is a non-dimensional ion charge 
density, and within the sheath n_#n , whereas within the 

A A A A i G 

ambipolar region n =n and n <n and in the free stream 

^ e i e Saha 

plasma n =n =ft . In fact, for purposes of this work, the 
e i Saha 

edge of the sheath will be where n approaches n within 

® i A 

1 %, and the edge of the ambipolar region will be where n 
^ e 

approaches n within 

Saha 



3* Pois so n* s Equation 



Using the non-dimensionaliziag scheme on Poissuu’s 
equation. Eg. 3 becomes 



V 




n . 
1 



( 3 ) 



^ • Energy Eq ua ti on 



E = 



Applying the non-dimensional parameters of Eg. 7 to 
-Vt gives the following result: 



(kT ) e 

E V(f) 



( 9 ) 



In the current equation, account must be taken of 

the effect of varying electron temperature. Since J >>J 

e i 

the current is essentially an electron current. Neglecting 
magnetic field effects and using Einstein's relation, the 
current is written as 
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( 10 ) 



kT T 

= -y^e [n.Vcj) - — — Vn ] 

e e ^ e T e 

o 



With the introduction of a dimensionless temperature. 



0 = T /T , the non-dimensionalized current equation becomes 
® ° y^(kTo)^ e/ 



J = 
e 



o 



8 



( 11 ) 



where 



/N /S /N 



D = -n^V(j) + eVn^ 



( 12 ) 



Equation 11 will contain the term (1+/3^) in the denominator 
of the coefficient when magnetic field effects are included. 
The dot product of Egs. 9 and 11 yields 



. SI 

n — 



(13) 



Non-dimensionalizing n and letting d ^ in PJq. 6, 

e s ^ rrig 

and introducing Eg. 13 into Eg. 6 results in 

p (kT )3 e 2 ^ ^ 

- _e 2^-2- j.y* = <’“> 

3ae 



As the electron temperature increases, the collision 
frequency with neutral particles increases affecting the 
terms fj. and Q . Since electron-neutral collisions are by 

far the most numerous Q can be approximated by Q= V 5 

en n 



or 



“ = "nOen«n<V>"n> 



8kT 



irm 



o gl/2 



(15) 



In the range of interest of the parameters, the collision 
cross-section changes little with temperature so Q depends 
1/2 

on e 
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A modification to Compton's eguation[27] for small 
values of E/P gives an electron mobility of the form 

- 1/2 






60 





/8kTo' 


nie ' 


\Jjfm ^en^n 



(16) 



-1/2 

so that jj gees as 6 . Combining Egs. 15, and 16 and 

e 

using the perfect gas law results in the form of the energy 
equation used in this work. 



Yj-V(}) = -n^(e^-e) 



(17) 



where 



Y 



0.13(ky^ e/ 

P^Q 6 
en e n 



(18) 



The value of the coefficient y appears to dictate 
the degree of Joule heating in the system. This can best be 
seen by expanding Eg. 17 and expressing 6 iraplicitely in the 
form 



6 



(y/e) Vd)- V(^ + 1 

(Y/0) vJ-Vn^/n^ + 1 



(19) 



Fcr physical reasons it is expected that 6>1. It 
can be seen from Eg. 19 that for very small values of y, 
0 — ^■l, which is the constant temperature case. Similarly, 
for y which is very large, an upper limit for 0 is found. 



e 



n^V$ • V4> 



for Y “ 



( 20 ) 



V$ • Vn^ 

It is interesting to note that for very large y, 6 becomes 
independent of y, which means that it is independent of the 
chemistry of the plasma. 



Now that y has been described, it will be useful to 
see what various HHD plasma compositions it represents. 
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Assuming atmospheric pressure, and the perfect gas law, Zg. 
18 can be shown to be 

Y = 10"^® X -|2_i (21) 

O Q ^ 
n en 

in MKS units where Z is the atomic weight of the principal 
neutral species, or a mass average of the neutral 

constituents . 

For combustion type gases where the principal 
constituent might be a typical value of y can be found. 

From Ref. 26, pages 139 and 148 for a temperature of lOOQOK, 

Q = 1.5x10-i«m2, d = 2x103, y is found to be 

en n 

98.3. This type of gas has y values which are low, 
characteristic of plasmas with high collision 
cross-sections. 

In many inert gas plasmas nitrogen is used to help 

the plasma reach equilibrium more quickly. The same 

reference gives typical values of Q = 5 x 10 ~ 2 Or [.2 and 6 = 

es n 

8. Then Y = 1.4x105. This is already at the upper limit 

described by Eg. 20. Even higher values are obtained with 
argon and cesium seed to the order of 10^ because of the 
Ramsauer effect[Ref. 27, p. 31]. 

5 • S_pe ci^ Eg ua t ions with Magnetic Fi e1^ Ef f ec t s an d 
Negligible Joule Heating 

When Joule heating is small and collisions between 
electrons and neutral particles are sufficient for efficient 
energy exchange, electron temperatures will not differ 
significantly from the neutrals. Furthermore, ion 
temperatures will be essentially the same as the neutrals 
regardless of Joule heating because of their large mass. 
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That is, the system of equations is stable even when T >T 

e o 

because it has been assumed that T =T = a constant which is 

o i 

unaffected hy Joule heating. The following expression 
replaces the energy equation for this case: 



T 



e 



T. 

1 



T 



O 



(22) 



Now Eg. 1 is solved for J and J in a two-dimensional 

X y 

Cartesian system with x and y coordinates. 

9n 



= (1+6^) ^ 



H * I’ 






(23) 



Jy = (l+B^)"^ 



- ^s^ 9y- 



+ g (-y^n e ± D 
s s s 9x 



1 

s® 9x ^ ^ 



(24) 



Substituting Eqs. 23 and 24 into Eq. 3 and using Egs. 5 and 
22 to help simplify results in 



-n V^4> 
s ^ 



V(})-Vn ± (kT /e)V n 

^ s ' o s 



+ 



9n^ 

^s 



9i _ ^^s 
9y 



9(|) 



9y 9y 



) 



= 0 
(25) 



Equation 25 represents two equations, one for electrons 
(+sign) and one for ions (-sign) . Now / 3 , - 0 because of the 

high ion inertia. Hence the two equations representing 
species continuity and Ohm's law become 

-n^V% - y4,-Vn^ + (kVe)v\ + g - -^ ||) - 0 (26) 

-n^V^4) - V(|)-Vn^ - (kT^e) = 0 (27) 

Non-dimensionalizi ng and substituting Eg. 8 into 26 and 27 
results in 
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ne(%-n^) - ^$*Vrig + + Bijj = 0 



/S /s /s 

"i<"i-''e> - 



sAy ^ A\ 



^n± - 9 n. = 0 



(28) 

(29) 



where 



/S /\ /s /s 

3n 3<t) 3n 3(J) 

. _ , e __ _ e , 

' /N /S As. As. ^ 

3x 3y 3y 3x 
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6. Spe cie s Equ ations with Non- C on st an t Elect ron 
Teinperat a re 

For the case of significant Joule heating Egs, 23 

and 29 are net descriptive of the species concentration 

since T #T everywhere in the field. It can be shown that 
e o 

for this case the species equations become 



/N /S 
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28 and 29. 



C. BOUNDARY CONDITIONS 

Many boundary conditions can be hypothesized for wall, 
electrode, free stream, and upstream/downstream locations. 
Free stream conditions are actually those at the edge of the 
ambipolar region. They are taken to be zero potential and 
zero space-charge, and a Saha-predicted total charge. The 
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electrode node potential is a fixed value. Other wall 
boundary conditions are not so easily determined since they 
are frequently vague and not well identified in the 
literature. Also, since this model is new, the boundary 
conditions used must be compatible with the system that is 
being represented. 

Successful runs have been generated using the following 

conditions with T = T ; 

e o 

1) Free stream as mentioned above, 

2) Periodic electrode nodes allowing periodic conditions 
at upstream and downstream interfaces (Pig. 1)/ 

3) Line electrode nodes allowing current constrictions 
(the necessary condition for two dimensions [ 2 1 ]) . The 
potential is fixed and ion and electron concentrations are 
zero at the node, as dictated by a catalytic surface, 

4) The inactive portion of the electrode wall is treated 
as an "insulated wall". Current perpendicular to this 
wall is zero, and the wall eliminates space charge and 
stagnates charge motion. 

In reality a small sheath will form along the insulated 
wall because of the difference between electron and ion 
mobilites. The analysis of Appendix C shows that the 
thickness of the sheath goes approximately as the square 
root of the potential. Since the floating potential will 
generally be much less than the electrode node potential 
drop it is expected that the insulated wall sheath will be 
very much smaller in size than the electrode sheath. 
Consequently, the hypothetical boundary conditions at the 
wall specified above are representative. 

A closer look at the electrode node boundary conditions 
is needed. Blue and Ingold[28] point out that though most 
authors use zero charge density at the electrode surface. 
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some cases 



for the 



require the modification n = 2J/(ec) 

s 

electrode boundary conditions. Appendix C gives 

characteristic values for J and c of 7x10'* amps/in^ and lOs 

m/sec, respectively, at typical temperatures. This results 

in n = 10*9 particles/m3 which is of the same order as the 
s 

charge densities which will be encountered later in this 
work. The difficulty in the application of this boundary 
condition comes from the fact that a true "point" electrode 
node would have an infinite current density. As will be 
seen later the numerical procedure uses a finite-difference 
scheme and the electrode is not a true "point". The current 
density at that location becomes artificially dependent upon 
the mesh-size. To avoid this problem, either the charge 
density must be set to zero, or some new non-zero restraint 
must be found. The use of periodic boundary conditions 
seems indicated since solutions of Laplace's equation 
(Poisson's equation in the absence of a space charge) tend 
to be periodic in two dimensions. As will be seen, several 
conditions were tried and the results are discussed in later 
sections. 

Other boundary conditions, such as a catalytic insulated 
wall were attempted with no reasonable results. This is 
discussed further in the next section. 



Figure 2 is a repres 


entation 


of the 


domain in 


the 


vicinity of the electrode. 


and 


some 


boundary 


conditions 


that 


are typical of those used 


in 


this 


work . 


The domain 


is 


represented numerically 


by 


a two-dimensional array 


of 


equally spaced points wh 


ich 


are 


operated 


on by finite 



differencing the controlling equations. 

Boundary conditions for $, n , and n remain the same as 

e i 
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above for the Joule heating case. Temperature boundary 
conditions are straightforward. End conditions are again 
periodic and free stream requires that 9=1. Since both 
the electrode node and walls consist of dense materials with 
which electrons collide, it is assumed that the electrons 
give up their surplus energy readily to these surfaces. 
Therefore, 9=1 along that boundary. 
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Figure 2. Arrangement of electrodes and boundary conditions 
for computational array. 



38 



% 



i 



III. hsthod of solution 



A. INTRODUCTION 

The governing equations (namely Egs. 8, 28, and 29) are 

solved numerically for the parameters n , and n 

e i 

throughout the domain as shown in Fig. 2. Later when the 
electron energy equation (Eg. 17) is included, a fourth 
parameter 0 is included in the field, with Egs. 31 and 32 
representing the species equations. The method of solution 
will now be discussed together with the application of the 
boundary conditions. The system of partial, non-linear, 
second-order differential equations is basically elliptic 
and constitutes a boundary value problem. The partial 
differential equations are replaced by finite difference 
equations; the technique for solution is discussed below. 
The non-linear nature of the controlling equations requires 
rather sophisticated computer techniques, and the boundary 
conditions present special problems. 

There are a number of possible approaches to the 
numerical analysis of the equations and several were tried 
before a successful approach evolved. In order to better 
understand the nature of these complicated equations, this 
Section also discusses the unsuccessful approaches and 
probable reasons for failure. 

B. NON-LINEARITY CONSIDERATIONS 

The two species equations contain the non-linear terms 
which present special programming problems. One-step 
techniques for solving this system of three simultaneous, 
second order, non-linear, partial differential equrtions are 
non-existent. One solution procedure in such cases is to 
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include all non-linear terms on the "right hand side," that 
is, external to the coefficient matrix, and hope that these 
non-linear terms change slowly enough with each iteration to 
render a convergent process. This procedure is known as trie 
Jacobi method[ 29 ]. Looking at equations 28 and 29 for the 
constant termperature case, or 31 and 32 for the Joule 
heating case, it can be seen that only one terra in each, 

r-7 2 ^2 

either V n or V n / is linear, and that there is no 
e i 

prescribed "load" on the right hand side. As will be seen 
in the expansion of these equations later, this leaves a 
dozen or more non-linear terms for the right hand side. 
Experience shows that the Jacobi method proves to be 
unstable . 

As a consequence of the failure of the Jacobi method, a 
quasi-Jacobi method is used. When the product of two 
variables is encountered, one variable is treated as a 
constant coefficient for each iteration. This means that 
the non-linear terms are retained in the coefficie'nt matrix. 
The "constant" coefficients are updated after every 
iteration, thus changing the coefficient matrix- 
Successively higher voltages were calculated in increments 
allowing small changes to the system. 

Had this quasi-Jacobi method not succeeded, the next 
step would gave been to apply a Newton-Eaphsen cechnigue[ 29 ] 
to the equations. Because of the increase in the number of 
terms and the consequent increase in computer requirements, 
the x'lewton-Raphsen method was reserved as a last resort. 

C. GRID SIZE AND COMPUTATIONAL SEQUENCE 

The number of points required to define a domain 
sufficiently large to show the sheath and ambipolar regions 
is critical (See the section on numerical results) . For 
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example, if a two-dimensional field consisting of a square 
array of 31 x 31 points is desired to be solved 
simultaneously, then it can be shown that a finite 
difference scheme will require in excess of 66 million bytes 
of storage for the coefficient matrix alone. Special 
storage techniques which take into account the banded 
structure of the coefficient matrix reduces the storage 
requirements to 12 million bytes. Even this is about 12 
times the core capacity of the IBM 360/67 at the Naval 
Postgraduate School. Obviously, to have a sufficiently 
large array and to keep within the core capacity the 
solution method must consider portions of the array at one 
time. 

This problem is met by using the line iterative method, 
i.e. by solving one line of the array at a time, sweeping 
back and forth with successive field iterations. By doing 
so, the size of the coefficient matrix is reduced to a small 
fraction of the full matrix. Because the equations are 
elliptic, and therefore each point affects every other point 
in the field, it is desirable to reduce the computer time 
for "information" to move through the field. This is 
accomplished by rotating the field 90° after every second 
field iteration such that the line sweeps go forward, then 
backward, then up and then down. This approach takes 
considerably less execution time for a converged solution 
than the usual procedure without alternating directions and 
rotation . 

The final grid size chosen for this line iterative 
method is 51 x 51. This decision was based upon the time 
required to reach a converged solution. The time for each 
run was of the order of four or five hours. The coefficient 
matrix requires only 51 x 3 x 10 x 8 = 12,240 bytes while 
the total space needed is 150,000 bytes for the constant 
temperature case, and 170,000 bytes for Joule hearing. This 
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storage includes program space, variable storage and 
solution matrices for double precision arithmetic. 



To further increase the speed of convergence a method of 
successive over-relaxation is incorporated similar to that 
outlined in Ref. [30]. At the end of each line iteration 
the parameters are subjected to the equation 



A =(l+w)A -wA,, 

corr new old 



(33) 



where A is the corrected 

corr 

the most recent calculation, 
w is the over-relaxation para 



value of the parameter, A is 

new 

A is the previous A ,and 
old corr 

meter. 



D. CONTROLLING EQUATIONS 



The finite difference representation of all terms was 
expressed to the order of when practical, where h is the 
grid spacing of the square 51 x 51 array. Ketter and 
Prowell (Ref. 29, p. 226-228) also describe the procedure 
for setting up finite difference equations. Second 
derivatives in two-dimensions are given by 






i,j ^^+1, j ^i-1, j ^i, j + 1 ^i/ j-1 



-4A. .)/h^ (3^*) 

1/3 



First derivatives have the form 
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(A. 



1+1 /j 



A. , J/(2h) 

i XfJ 



(35) 



Where a first derivative is required at a boundary it is 
written as either a forward or backward difference. 



9A 

9x. . 
1/3 



= (-3A. - A.^2_.)/(2h) 



(36) 



Transforming Eg. 8 into finite differences gives 

/s /V ^ 

+1+1,3 ^1-1/3 ^1/3+1 ^1/3-1 1/3 



2^ 
+ h n. 



1 . . 
1/3 



( 37 ) 



7'' 
- h^n 



= 0 



1/3 
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This is the simplest of the three controlling equations. 



The species equations are considerably more complicated. 
The finite difference form of Eq. 28 is transformed here to 
illustrate the use of variables as "constant" coefficients 
in the non-linear terms. 







-n 

1 ''i+l,j ’ 


^ <"e. 


/\ 

+ n 

e . . 


- 4n^ )/h 



e . 
1 

2 



) - n + n 

1 ®i,j+l ®i,j-l 






' "e. .‘"e. .-"1. •> + - * 1 - 1 , j> .'"e. , . 

ifl 1/D i/D 1+1 /D 1-1 /D 



+ 6(n -n )]/4h - (n +n ) /h (38) 

i,j+l ®i,j-l ^i+l/j ^i-l/j 



This equation corresponds to a line in the y direction. The 

bracketed terms {} are treated as constant coefficients but 

are updated with each field iteration. The finite 

difference version of Eg, 29 resembles Eg. 38, but is 

considerably simpler since =0. 

i 



When the electron temperature is uniform and equal to 

the neutral temperature, no energy equation per se is 

required since the condition T = T is incorporated in the 

e o 

two species equations. However, for the non-equilibrium 
case, the energy equation is solved separately from the 
other three. That is, Eqs. 8, 31, and 32 are solved 

simultaneously for assumed values of 9. With each iteration 
6 is updated by solving Eq. 17. 

E. BOUNDARY CONDITIONS 



Free stream boundary conditions are the simplest. For 

each case the free stream values are set at ^ = 0, n = n = 

e i 

n , and 0 = 1 . 

Saha 



Upstream/downstream conditions are chosen to be 
periodic. Ihey require that not only the values of the 
unknowns be the same at periodic points, but also that their 
derivatives be the same. For example, if the far upstream 
station is labeled 'M'* and the far downstream station is "a" 
the two equations, then 



^ . = A, . and , 

n,3 1,J 



A, . = — A . + A T . 

Ifl 2,3 n ,3 n-l,D 



(39) 

(40) 



would apply for all parameters $, n , n , and 0. 

e i 



The insulated wall is hypothesized to be neutral and to 

equilibrate the charges, therefore, the conditions are 

numerically set to n = n = n . The wall condition is 

e i Saha 

taken from the restriction that perpendicular current is 
zero. From a non-dimensional form of Zq. 24 for zero 
current in the y-direction the boundary condition at the 
wall is 

3n 



-n. 



3$ 


✓s 

8n 

0 


3$ 


/s 


/N 


- 3n -^ + 3 


3y 


3y 


® 3x 



= 0 



8x 



(41) 



The x-derivative of n is zero along the wall and n = 



A 

n 



Saha 



So 



n 



saha 



8(|) 

3^ 



3n, 



3(j) 



^ ^^saha 
3Y 3x 



(42) 



The finite difference expression for Eg. 42 is 
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( 43 ) 



The boundary conditions at the electrode node consist of 



with catalytic electrodes. 



F. EVOLUTION OF THE MODEL AND PRIOR ATTEMPTS 

The factors affecting the choice of the calculation 



conditions, array size, electrode placement, choice of 
coordinate system, and the form of the controlling eguatioas 
are but a few of the things that had to be considered in 
modeling the problem. The final successful technique 
evolved by trial and error. The following paragraphs give a 
brief account of the evolution including successes and 
failures . 

At the outset a two-dimensional, flat plate, continuous 
electrode model was proposed. Because of the small 
thickness of the sheath relative to the electrode size, most 
areas of the sheath would not experience end effects, and 
would be effectively one-dimensional. Further investigation 
led to the proof of the non-existence of a one-dimensional 
solution presented in Appendix B. 

It then became obvious that the key to this problem was 
in two- or three-dimensional current constriction since the 
one-dimensional Cartesian geometry offers no natural 
geometric means to decrease current density away from tne 
electrode. The electrode is then represented by a series of 
nodes as shown in Fig. 1. An examination of the literature 



the three conditions , n =0, and n 

o i 



e 



0 consistent 



procedure and computational model are many. Boundary 
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showed that this mode of current constriction was known and 
even helps explain the electrode spots which may be present 
on the anode. 

At first the model consisted of both the anode aiid 
cathode on opposing sides of a field of computational 
points, but it was realized that a field of this size could 
not possibly show all the details of the sheath. 
Conversely, a field capable of showing the required details 
could not contain both the electrodes. So the final model 
consisted of a series of line electrodes (for simplicity, 
the anode) along an "insulated" wall in a periodic field. 
The opposite boundary would then represent a free stream 
conditon. The other two boundaries are u pstream/downstreara 
boundaries. 

Once the model was chosen, the job of solving the 

equations was undertaken. Initially, the three controlling 

equations (the energy equation had not yet been considered) 

were solved separately for n , and n , respectively, in 

e i ' 

iterative fashion. The equations were Eqs. 3, 26 and 27 
with a Hall parameter of zero. The calculations diverged 
almost immediately because in the vicinity of the electrode 
many derivatives are large and values of the parameters 
changed too rapidly from the initial "guessed" values. 

It became apparent then that the three equations would 
need to be solved simultaneously. Section III.C explains 
some of the problems of such a technique in terms of overall 
storage requirements. In the light of this, the line 
iterative method was chosen. 

Armed with a simple set of variables and the 
line- iterative method a solution was attempxed. Convergence 
was achieved, but very slowly, since it took many iterations 
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array was 



for "information” to cross the array. So the 
rotated 90° at intervals causing the effective sweep to be 
back and forth, then up and down. 
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A variable mesh size was also attempted 
success, but since it was determined that the 
region was of the same order of magnitude as the s 
was rejected as an unnecessary complication. 

Various boundary conditions were tried, but mo 
did not work or gave physically meaningless resul 
most difficult boundary to model is the insula 
From the final boundary conditions chosen the s 
equations describes a certain physical case. 

Once the present technique was successful for 
voltages, it proved its versatility by accep 
addition of a magnetic field and the energy 
routinely. While not being the ultimate in 
techniques for solving this system of equations, th 
scheme has rendered useful results and promis 
further productive in other cases. 
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IV. RESULTS 



A. PROCEDURE AND CONVERGENCE CRITERIA 

The simplest case, which included no magnetic field 
effects and negligible Joule heating began with the 
following assumed values: 

1. zero potential everywhere except at the electrode, 

2. selected non-dimensional potential at the electrode, 

3. zero charge density everywhere except at the 

boundaries, and 

4. Saha charge density at the boundaries. 

Using these assumed values, iterations for improved values 
proceeded . 

Periodically convergence was checked after many field 
iterations by comparing the left hand side of Eg. 8 with the 
right hand side at various locations in tne field. This 
convergence criterion may appear somewhat arbitrary, but was 
chosen when it was found that Eg. 8 was more difficult to 
satisfy than Egs. 28 and 29 and therefore represented 
stricter convergence criterion. Additionally, Eg. 6 

contains all three unknowns providing a check of a3 1 
parameters. 

The error was calculated from the finite difference 
representation of Eg. 37. 



As the calculations converged the error decreased with 
subseguent iterations, but the rate of convergence slowed as 
the number of iterations became large. Eventually a point 



6 = 1 



([) . . _ + (}) . . , + 
^ 1 , 3+1 ^ 1 , 3-1 

V.2 




(44) 
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was reached where further iterations were not cost eff 
as the additional improvement in the solution re 
excessive computer time. A second convergence criteri 
used in addition to the the first. Certain potential 
at critical points in the array were plotted witn succ 
field iterations. If these values approached 
assymptotic value, then the procedure was cons 
convergent. 



ective 
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values 
essive 
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idered 



As discussed earlier, a successive over-relaxation 
technique was used to increase the rate of convergence. It 
was found that over-relaxation parameters exceeding 1.1 
caused the field to diverge for all but very low voltages. 
Some improvement of convergence rate was noted for w = 1.05 
so this value was used for most runs. When attempting to 
achieve solutions for higher voltages or other conditions 
which might prove divergent, such as large Hall parameters 
or Joule heating parameters, a value of w = 0 was used in 
order to put the least possible "stress" on the system of 
equations. 



B. NEGLIGIBLE JOULE HEATING AND NO MAGNETIC FIELD EFFECTS 



Convergence for this set of runs was slow. In 
achieve an error of less than 20% (|6|<.2) 

iterations were required, even with a near 
over-relaxation parameter of 1.05. This equates 
four hours of computer on the IBM 360/67. An a 
hour (250 iterations) reduced the error to 16 
results shown in Figs. 3-13 were completed with 
Further improvement would have required excessive 
time. 
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Figures 3, 6, 9, and 12 are plots for the case 



I = 



5.80. 
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Figure 3. Potential^ and Charge Density Figure 4. Potential and Charge Density 

Profile, 5)=5.80, 7=0,/9=0, n=10“-^, Profile, 0=17.4,7 =0,/3 =0, n=10" 

array size=500L, 51x51 grid. array size=500L, 51x51 grid. 
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Figure 7- Potential Contour Plot, (^17.^, Figure 8. Potential Conto^ Plot, $=29.0 

y=0,/3=0, n=10"3, array size = y=0,/3=®> ^=1*^ > array size = 

500L, 51x51 grid. 500L, 51x51 grid. 
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Figure 11, Electron Density Contour Plot, Figure 12. Ion Density Contour Plot, ^5.80, 

$=29.0, =0, n=10"3, array y=0>/3=0, n=10~ , array size=500L 

size=500L, 51x51 grid. 51x51 grid. 
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Figure 15. Potential and Charge Density • Figure l6. Potential Contour Plot, $=29.0 

Profile, J)=29.0, y=0, ^=0, h=10"^, y=0,/3=0, n=10“3^ array size= 

array size= lOOOL, 51x51 grid. lOOOL, 51x51 grid. 
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Figure 17. JElectron Density Contour Plot, (f)=29.0. Figure l8. Ion Density Contour Plot, ^29.0, 

y=0,(3==0, n=10-3, array si7,e=> tOOO L, y=0,/?=0, n=10"'^, array size=|000L, 

51x51 grid. 51x51 grid. 
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This value represents one volt at T =200Qok. Figures 4, 7, 

o 

10, and 13 are for $ = 17.41 and Figs. 5, 8, 11, and 14 are 

for $ = 29.02, representing 3 and 5 volts, respectively, at 

T = 20000K. 
o 



The overall size of the array for each of these runs is 
500 times the characteristic length defined in Eg. 7. At 
2000OK this means that the dimensions of the sgiiare field 
shown in the figures are 5.25x10-sm on a side. The 
non-dimensional charge density as determined by the boundary 
condition is 10~3. This corresponds to 8.5x10i7 
particles/m3 , which is typical of an equilibrium charge 
density for seeded plasmas at this temperature. 



1 Di str ibuti on s 

Figures 3-5 are all one-dimensional profiles of 
potential and charge density distributions along a line 
perpendicular to the electrode wall. The aoscissa shows 
distance from the electrode in terias of the characteristic 
length. Figures 6-8 are two-dimensional contour plots of 
the potentials for the three voltages. The electrode is 
located at the bottom center where the contours concentrate. 
Except for the absolute magnitudes, there is litrle change 
noted in the potential plots for the three voltages. As 
might be expected, the steepest potential slope occurs near 
the electrode where the greatest electric field exists. At 
the free boundary the slope of the potential is nearly 
constant indicating a constant electric field. A constant 
electric field is a property of two opposing infinite 
flat-plate electrodes. The plasma spreads the effects of a 
point or line electrode so that the free stream behaves as 
if the electrodes were infinite plates. 
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Ch arged Pa r t id e Di st ributions 

Figures 9-11 demonstrate considerable change for the 
electron densities with increasing voltage as do Figs. 12-14 
for ion densities. 

An examination of the charge density plots in Figs. 

3f 4, and 5 reveal several interesting facts. First, siiice 

the edge cf the sheath is the line of points where the 

charge densities are equal, the figures put the sheath 

length at 35%, 90%, and 1005 of the total length of the 

field. In the following table, the characteristic sheath 

length as calculated from Eg. C.7 for T = 200QOK is shown 

o 

next to the computer generated sheath length from the 
extimate cf the convergence of the charge density lines in 
the figures: 



$ 4> cliaracterist ic computer 

o o 









sheath length 


gen era ted 


5.80 


1.0 


volts 


-5 

.800x10 m 


-5 

1.84x10 


17.41 


3.0 


volts 


-5 

1.38x10 m 


-5 

4.72x10 


29.02 


5.0 


volts 


-5 

1.79x10 m 


-5 

5.25x10 



These offer order-of-magni tude comparisons at best, but the 
predictions of the characteristic sheath length with Eg. C.7 
are useful in estimating the size of the field for the 
computer program. 

Notice that the charge density slopes on the right 
of Figs. ^ and 5 are significantly greater than zero. This 
means that the field was not sufficiently large to encompass 
the entire ambipolar region. Since the boundary conditions 
at the free stream are imposed as con.stants they are 
'•forced" conditions and probably result in error for these 
cases. One result of this is the erratic behavior of the 
charge density close to the electrode. As will be seen 
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later, this behavior has little effect on the current 
distribution as long as it occurs in a confined region. 

It is evident that the observation of 100% of the 
field for the sheath length in Fig. 5 is probably low since 
the convergence of the two charge lines occurs artificially 
at the boundary. To prove this, and the other observations 
made above, a run was made for the $ = 29.02 case doubling 
the field length to 1000 times the characteristic length, or 
10.5x10~5 in. Figure 15 shows the dramatic result. The 
erratic behavior of the charge density lines near the 
electrode is gone and the slope is essentially 2 ero near the 
free boundary. The sheath length is measured at about 
3.94x10-5 m, 2.3 times the characteristic predicted value. 
It is interesting to note that this is about the same ratio 
as the measured/predicted ratio of the $ = 5.80 case. 
Figures 16-18 show the two-dimensional plots for this 
$ = 29.02 case. 



C. NON-CCIISTANT ELECTRON TEMPERATURE 

With the addition of the energy equation the effect of 
Joule heating on the system was studied. The value Y = 100 
was introduced first to simulate combustion gases. The 
results are shown in Figs. 19-22 for the cases $ = 5.80 and 
29.02. In spite of the earlier discussion showing that the 
size of the field is inadequate for the higher potential 
case, it is used here to avoid the complication of changing 
the electrode spacing when comparing different potentials. 

The most obvious change caused by the heating cf the 
electrons is the changed slope of the potential indicating 
less shielding for the Joule heating case. This is further 
shown in Figs. 23 and 24 which represent $ = 5. 80 and 
y= 105. The charge density profiles are relatively 
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later, this behavior has little effect on the current 
distribution as long as it occurs in a confined region. 

It is evident that the observation of 1005J of the 
field for the sheath length in Fig. 5 is probably low since 
the convergence of the two charge lines occurs artificially 
at the boundary. To prove this, and the other observations 
made above, a run was made for the $ = 29.02 case doubling 
the field length to 1000 times the characteristic length, or 
10.5x10~5 m , Figure 15 shows the dramatic result. The 
erratic behavior of the charge density lines near the 
electrode is gone and the slope is essentially zero near the 
free boundary. The sheath length is measured at about 
3.94x10-5 m, 2.3 times the character istic predicted value. 
It is interesting to note that this is about the same ratio 
as the ireasured/predicted ratio of the $ = 5.80 case. 
Figures 16-18 show the two-dimensional plots for this 
$ = 29.02 case. 



C. NON-CONSTANT ELECTRON TEMPERATURE 

With the addition of the energy eguation the effect of 
Joule heating on the system was studied. The value ~Y = 100 
was introduced first to simulate combustion gases. The 
results are shown in Figs. 19-22 for the cases $ = 5.80 and 
29.02. In spite of the earlier discussion showing that the 
size of the field is inadequate for the higher potential 
case, it is used here to avoid the complication of changing 
the electrode spacing when comparing different potentials. 

The most obvious change caused by the heating of the 
electrons is the changed slope of the potential indicating 
less shielding for the Joule heating case. This is further 
shown in Figs. 23 and 24 which represent $ = 5.80 and 
y= 103. The charge density profiles are relatively 
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unaffected by the addition of high temperature electrons, in 
fact there is no observable change in the sheath length. 

Figure 25 is the electron temperature profile (T /T vs. 

e o 

distance from electrode) using the same cut as the previous 
one-dimensional plots. This plot shows that the majority of 
the Joule heating occurs very close to the electrode where 
the largest electric fields exist. There is no noticible 
electron temperature rise in the arabipolar and free stream 
regions. This graph is typical of the temperature profiles 
for the various cases and differs only in the height of the 
peak near the electrode. This peak appears to be a function 
of the potential as well as the Joule heating parameter. 

From the numerical standpoint, the addition of Joule 
heating offers a benefit not easily foreseen. Convergence 
of the system with the addition of the energy equation vas 
accelerated such that errors of less than 2% ware achieved 
in less than 250 iterations. The reason for this 

improvement is not yet clear but its advantages are obvious. 

The inclusion of Joule heating is a better approximation for 
most generators, and provides the bonus of reducing the tirae 
of numerical convergence. 

d: effect of equilibrium density 

An interesting case is that of increasing the 
equilibrium charge density to much higher values by changing 
the boundary value of the electron and ion charge densities. 
Presented here are the results of setting ft = 1 at the 

boundary. At a temperature of 2000°K this corresponds to a 
charge density of lO^i particles/m^ . Though not a common 

case, it serves as an indicator of some of the numerical 
properties of the method, and a good comparison for the 
sheath length prediction. 
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Figure 23. Potential^and Charge Density Figure 24. Potential Contour Plot, ^5.80, 

Profile, (J)=5.80, y=1000, ^=0, n=10"3^ "yalOCO /3=0, n=10“^, array size 

array size=500L, 51x51 grid. 500L, 51x51 grid. 
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Figure 25. 



Electron Temperature Profile, $=5.80, y=100, 
/3=0, n=10“'^, array size=500L, 51x51 grid. 
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At T = 200QOK and for $ = 5.80 and n = 1 at the 

o 

boundary, the predicted sheath length from Sg. C.7 is 

2.5x10~7 m, At this same temperature the characteristic 
length as defined in Eg. 7 is L = 1.0x10-^ m. Figures 26 
and 27 show the case where the field length is 25L or 

2.62x10“* m . As expected the sheath length is about 10?, of 

the total field. 



The dotted portion of Fig. 27 is drawn to agree with 
previous arguments on the electrode boundary conditions. 
Actually, some of these higher density runs were made using 
a different boundary condition, namely one that attempted to 
impose a consistent total current through the field. 
Though erroneous, it was found that this particular boundary 
condition has little effect on the rest of the field, and 
that conditions were actually determined by the 
non-conducting wall surface. The "kinks" seen close to the 
electrode in several of the one-dimensional plots appear to 
be a resolution feature of the numerics which tends to 
override the effects of the charge density boundary 
condition at the electrode. 



E. EFFECT OF THE NUMBER OF ELEMENTS IN THE ARRAY AND 
ELECTRODE POSITIONING 



Earlier trials used a 41x41 array instead of the 51x51 
array used later. In general the smaller arrays did not 
encompass the entire ambipolar region as can be seen in 
Figs. 28-30 which are for n = 1.0, cj = 11.6, and array 
size = 20L, 25L, and 35L, respectively, in 41x41 arrays. 
They seem to lack sufficient room to satisfy the equations 
even after varying the mesh size. 



Figures 28-30 also show the changes that occur when the 
mesh size is changed. The sheath length appears to be 
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Figure 26 . Potential and Charge Density Figure 27 . Potential and Charge Density 

Profile, $=5.80, y= 0 , / 3 = 0 , n=l. Profile, ^= 17 . 4 , ^= 0 , ^= 0 , n^ 

array size= 25 L, 51 x 51 grid. array size= 25 L, 51 x 51 grid. 



consistent with the scaling, but the ambipolar slope does 
not appear to change, indicating the inadequacy of 
mesh-sizing to control the size of the field when the number 
of elements is too low. 

For a qualitative look at a multiple electrode plot. 
Figs. 31-33 are provided. This scheme allows a current 
density comparison of this run with one of lower potential 
and the necessary symmetry in the absence of a magnetic 
field. Current profiles will be considered later. 

F. HAGNSTIC FIELD EFFECTS 

The introduction of a magnetic field induces the most 
visible, although not unexpected, effects of all the 
phenomena discussed so far. As might be anticipated, the 
symmetry of the field is lost as shown in Figs. 34-43. 
Notice however, that while the potential shows significant 
distortion, the charge densities do not, at least for the 
case of $ = 5.80. 

The angle of incidence of the potential contours with 

the insulated wall is approximately the complement of the 

Hall angle (arctangent ) as can be seen in Fig. 35 (450) 

and Fig 43 (270) . This is easily understood by solving Eg. 

24 with J =0 and considering only conduction effects. The 
Y 

result is 




This equation says that at points where J =0 and diffusion 

y 

is small, the slope of the potential lines is ~p. In the 
vicinity of the electrode where diffusion plays a more 
significant part the slope begins to deviate. 
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Figure 30. Potential and Charge Density Figure 31. Potential Contour Plot, $=29, ( 

Profile, $=11,6, ^=0, p =0, n=l, 'y=0.j3=0, n=l, array size=50L 

array size=35L, 4lx4l grid. 4lx4l grid. 
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Figure 32. Space Charge Plot (ni-iig) for multiple Figure 33. Total Charge Plot (n^+ng) for multiple 
electrodes, ^=2^.0,y-0, j2=0, n=l, electrodes, $= 29 . 0 , "^=0, /3=0, h=l, 

array size=50L, 4lx4l grid. array size=50L, 4lx4l grid. 
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Figure 3 ^. Potential^and Charge I^nsity „ Figure 35 . Potential Contour Plot, $=5,80 

Profile, <{)=5.80, y=0, p=l, n=10 y=0,^=l, n=10“ , array size= 

array size=500L, 51x51 grid. 500L, 51x51 grid. 
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Figure 36. Electron Density Contour Plot, Figure 37. Ion Density Contour Plot, 6=5.80 

$= 5 . 80 , y=0,^=l, n=10“'^, array y=0,/3=l, n=10“3^ eirray sise= 

size=500L, 51x51 grid. 500L, 51x51 grid. 
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0 2 50 500 250 0 250 

l/d ^ 

Figure 38. Potential^and Charge Density Figure 39. Potential Contour Plot, ^=29.0, 

Profile, $= 29 . 0 , y=0,/3=l, n=10“^, *7=0, /3=1, A=10“3, array size=500L, 

array si?,e= 500 L, 51x51 grid. 51x51 grid. 
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Figure UO. Electron Density Contour Plot, Figure 4l. Ion Density Contour Plot, ^29.0, 

$= 29 . 0 , y=0,/3=l, S= 10 " 3 , array 'y=0,|3=l, fi=>10“3, array si7.e=500L 

si?,e=500L, 51x51 grid. 51x51 grid. 
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FigTire k2. Potential and Charge Density Figure 43. Potential Contour Plot, ^5.80, 

Profile, $=5.80, y=0,;3=2, n=l, 7=0, /^=2, ft=l, array size=25L, 

array size=25L, 51x51 grid. 51x51 grid. 



As with the runs in the absence of a magnetic field the 
effects of the presence of a point electrode are dampened as 
the free boundary is approached. The magnetic field 
produces little or no distortion to symmetry outside the 
sheath. 



G. CURRENT PROFILES 

For each of the previous runs a complete two-dimensional 
current distribution can be made using the potential and 
charge density distributions. The current lines make an 
angle with the potential contour lines approximately equal 
to the Hall angle and would theoretically be exact except 
for the effects of diffusion. Figures 44 and 45 show the 
current streamlines in the field for ^ = 5.80 and P~ 0 and 
1, respectively. These were sketched from numerical data 
from the "Current Program" vfhich gives the current in vector 
form. The streamlines entering the field from the right are 
from the active site that lies periodically to the right of 
the one in the field. 

To determine the total current flowing through the 
system it is necessary to look some distance frcm the 
electrode. Theoretically/ for this model the density of the 
current will vary from some given value at the free stream 
to infinity at the electrode, however, as explained for the 
case of charge density a "smearing" takes place due to the 
numerical approximation. 

In terms of the amount of information available. Figs. 
46 through 48 are the most important in this work. They are 
the current-voltage diagrams for the charge densities 
= 10-3 and n = 1. The dimensionless current density was 
taken one or two stations from the free boundary (in order 
to avoid numerical errors at the boundary) and represents 
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Figure 44. Current Density Contour Plot, ^=5.80,y-0,/3=0, 
n=10~3j array size=500L> 51x51 grid. 
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Figure 45. Current Density Contour Plot, ^=5.80,y =0, /3=1, 
n=10“^, array size=500L, 51x51 grid. 
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the macroscopic current that would be observed in a 
two-dimensional model. The non-dimensionalization of tne 
current was consistent with Egs. 7 and 11 and contain the 
term ( 1 +^ 2 j. The conductivity was introduced to replace the 
mobility, and constant terms were collected to give the 
simplified coefficient shown in the figures. An 
interpretation of these results as well as examples of their 
use is given in Section V as conclusions. 
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Figure 48. Current- Voltage Diagram, n=l, 51^51 grid. 
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V. 



CONCLUSIONS AND RECOhU EN D ATIONS 



A. PHYSICAL CONCLUSIONS 

A successful model of the sheath evolved from the 
assumptions of steady state, frozen flow, and uniform 
temperature for neutral and ion particles. Qualitatively, 
the results of the analysis used in this work offer much 
insight into voltage drops attributable to electrodes in 
contact with an KHD plasma. Quantitatively, these results 
are useful if it is remembered that the model is 
two-dimensional and does not give the additional degree of 
freedom for current expansion afforded by three dimensions. 
A summary of the more basic conclusions which were drawn 
from the results is presented here. They will be explained 
individually later in this section. 



1 . 


The sheath 


can be 


self- gen era ted 


from a consis 


tent set 


of 


equations 


v/ith 


or without the 


use of 


a n 


energy 


equation. 












2. 


The electrode 


node boundary 


condition 


for 


charged 



particle density has little effect upon the field, while 
the insulated wall boundary condition has a profound 
effect . 

3. Current constrictions are necessary at the electrode 
to satisfy the system of equations. 

4. The resulting current- voltage diagram has a curvature 
consistent with theoretical predictions for space-charge 
in a one-dimensional flow of current and experimental 
evidence. 

5. The current density, for a given potential drop, 
varies inversely as the node spacing. 

6. The conductivity of the plasma within the sheath is 
critical in determining the sheath voltage drop, and 
methods to increase this conductivity will result in a 
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decrease of both sheath and boundary layer drops. 



1 • Sheath For ma t i on and Cha rge De nsi ty Pro fil es 

The results have shown that a sheath develops 
consistent with theory. The thickness of the sheath varies 
with the potential and is approximately that predicted from 
theory on Debye shielding. In each of the resulting plots 
presented in the previous section the size of the ambipolar 
region is essentially of the order of magnitude of the 
sheath length. No quantitative conclusions have been drawn 
as to the relative size of the ambipolar region other than 
to notice that higher voltages showed the tendency of the 
ambipolar region to increase in length relative to the 
sheath. It is very possible that at higher voltages the 
ambipolar region may grow to about the size of rhe boundary 
layer thickness. But since the potential drop in this 
region is low compared to that in the sheath it is not as 
significant in this analysis. The addition of Joule heating 
and a magnetic field did not appear to affect the sheath 
length. For the parameters used in this investigation, the 
difference between unlike charge densities (space charge) 
was about ' 10 % of the equilibrium density. 

Seme of the resulting one-dimensional profiles shown 
in the previous section exhibit erratic behavior in the 
first few stations past the electrode node. This is 
believed to oe caused by numerical convergence problems in 
the procedure arising from steep derivatives near the node 
and the extreme non-linearity of the equations. This 
appears to have little or no effect on the rest of the field 
since, as can be shown, the current distributions derived 
from these runs are consistent with those from smoother 
runs. Furthermore, charge densities other than zero were 
tested as a boundary condition at the electrode with no 
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change in the overall results. It seems then, that other 
conditions must be more dominant in the behavior cf the 
potential field. 

One very significant factor appears to be the choice 

for the insulated boundary condition. For example, a 

catalytic wall was attempted by setting n =fi =0, but no 

i e 

reasonable results were obtained. The biggest problem with 
the catalytic condition appeared to be the non-compa lability 
of the wall condition with the free stream boundary 
conditions, specifically the free stream slopes did not 
approach zero regardless of the mesh size or characteristic 
length. Density gradients occurring at a catalytic wall 
occur in a small region when compared with the sheath, and 
violate the assumption of a continuous fluid. Numerically 
this means that the charge densities must go from zero to 
equilibrium values within one grid space. Physically the 
majority of the field should not "see" this wall anyway. So 
the conditions chosen, namely that of an equilibrium wall, 
are consistent with the rest of the field and yield 
meaningful results. 

2 • C ur ren t Density Di stribu tion s 

The current- voltag e diagrams in Figs, hb through 48 
show the behavior of the system for various conditions. 
Figure 46 shows that the non-dimensional current density is 
inversely proportional to the electrode spacing since the 
current density for the $ = 29.02 case with the array 
size=1000L is approximately half that with the array 
size=500L. This is shown even more vividly in Fig. 47 where 
the introduction of multiple electrodes with double the step 
size (keeping the electrode spacing constant) resulted in a 
single current line. The significance of this dopendancy is 
that the voltage drop for a given current density becomes a 
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function of the electrode node density, and is thus a 
problen which has inputs from surface physics as well as HHD 
conditions. It must be remembered that for the s^juare 
arrays used in this work the electrode node spacinys are of 
the same crder as the sheath length. The minimum spacing 
must be of this order to avoid one-dimensional effects that 
would result from further "crowding" of the active sites. 
Spacing the nodes much greater than a sheath length makes 
computation more difficult since a much larger array would 
be needed to describe tne effects. Thus, this work models a 
minimum spacing resulting in a near maximum current density 
for a given potential drop. On an actual electrode, the 
spacing would be random and dependent on empirical factors. 

Figure 48 is the current-voltage diagram for the 
= 1 case using a 51x51 array. It was noticed that for a 
given potential, the current density for the P= 2 line was 
approximately twice that of the /3 = 1 line. Using this 
result a p> = 2 line was sketched in J'ig. 46, since no 
converged solutions were available for this line at 
n = 10-3. 



The introduction of the magnetic field moves the 
line to the right as shown in Fig. 46, and appears to cause 
a smaller potential drop for a given current. But recall 
that the current is already multiplied by the term (1+/3^) 
due to the effective reduction in conductivity, and this 
partially offsets the reduction in the potential drop. 

The introduction of the energy equation shows that 
Joule heating increases the current for a given voltage, but 
only to a slight degree. This is because the region that is 
affected by the Joule heating is small and very close to the 
electrode, well within the sheath. 

In the absence of a magnetic field the 
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current- voltage diagram tends to be concave upward, a fact 
that can be predicted with the following arguments: Cobine 
(Eef. 27, p. 129) shows that the high pressure space-charge 
current should go as V2/y3 in one dimension. If y is taken 
to be the sheath length in Eg. C.7 then it can be shown that 
V goes as J2. Figures 46 and 47 are nearly parabolic in 
appearance, and have curvatures consistent with this 
analysis. 

The zero voltage intercept produces an interesting 
problem since it has different signs for Figs. 46 and 47. 
It would be expected from Eg. 1 that for zero potential 
difference in the absence of a magnetic field, the electron 
current density would be positive because of the diffusion 
term and the positive charge gradient for rhe existing 
boundary conditions. However, recall that the validity of 
the catalytic electrode depends upon the fact that it is 
charged. In the absence of a charged electrode, the node 
should appear the same as the insulated wall, and no density 
gradient should exist. Therefore, no special significance 
can be placed on this intercept without further 

investigation of the boundary conditions. 

Reference [31] reports on the beneficial effect that 

increasing the seed fraction has on lowering the electrode 

voltage drops. From the behavior of the electron collision 

cross-sections with temperature, the optimum 

conductivity[ 32 ] is reached at seed fractions higher than in 

the core because of the lower total temperature at the 

electrodes. Heretofore, the only benefit would seem to be a 

higher boundary layer conductivity, but as shown in Figs. 

46-48, the local conductivity also affects the sheath drop. 

Hence, a higher value of O is expected to decrease the 

o 

voltage loss associated with both the sheath and ambipolar 
regions. 
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An example would be appropriate here to show the 
orders of magnitude of the quantities being dealt with. 
Suppose it was desired to find the potential drop across the 
sheath in a plasma where the temperature is I 8 OOOK, the 
current density is 0.2 amps-cm~ 2 ^ free-stream conductivity 
is 0.10 mhos/m, and the charge density is known to be of the 
order of 10^® particles/m® . Equation 7 shows the 
characteristic charge density to be 7 x 10^0 particles/m^. 
Thus ?ig. 46 would be the appropriate current-voltage 
diagram with n=10~3. The non-dimensional current density in 
the absence of a magnetic field is 1.50x10“® which gives a 
dimensionless potential drop of 32.0. This corresponds to 
4.9 volts at the given temperature. 



3 . Comparis on s wit h Experiment 

Argyropoulos et al[33] cite the use of a 

sophistocat ed computer solution for the boundary layer drop, 

comparing the results with an experiment on the AVCO-APL 

channel. Their results show the anode drop to be 101 volts 

due to the boundary layer. They do not consider sheath 

effects. Using these data, namely, T = 2000®K, 

wal 1 

J = 3.95x10^ amps/m 2 , /3 = 1.02 and a chemistry consisting of 
toluene and oxygen with cesium seed, the conductivity at the 
electrode is 1.14 mhos/m. Figure 46 shows a non-dimensional 
current density of 4.307x10“®. Using the /3 = 1 line this 
gives $ = 35.6 for an actual drop of 6.13 volts at the given 
temperature . 



Another comparison 
Teno[34]. For this experira 
constituents as the previ 
wall temperature of 1800®K 
parameters are J = 1.6x10 
gives j = 1.7x10“*, ^ = 121 



is from the work of 



ent 


they 


used the 


same 


ous 


one. 


The conductiv 


is 


0.30 


mhos/m 


and 
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sheath 


voltag 



Sonju and 
chemical 
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18.8 volts. 



Caution must be exercised when making these 
comparisons since several factors have not been considered 
yet. First, tendencies towards non-eguilibriu lo would 

increase the conductivity near the electrode resulting in a 
smaller potential drop. Secondly, nothing has been said 
about the node density of the experiments as compared to the 
assumptions made in the numerical solutions, namely that the 
distance between nodes is of the order of the sheath length. 
Should the actual node density differ from this it will 
affect the current-voltage characteristics. Thirdly, the 
use of a two-dimensional model rather than a 

three-dimensional one may tend to predict different voltage 
drops. The following example illustrates these limitations. 

An experiment by Kessler and Eustis[6] was conducted 
using ethanol in oxygen with 1% KOH. Nitrogen was 

introduced for cooling such that the N /O ratio was 0.5. 

2 2 

For this run T = 1685 ok, T = 270QOK, J=. 75x1 0 “^amps/ms , 
wall core 

( 3 = 1.5, and the electrode conductivity is 0.56 mhos/m. 
Figure 46 gives j = 3.697x10”'^, 4> = 490.2, and <{) = 71.1 

volts. This is a higher drop than the measured 45 volts 
that should account for both the sheath and boundary layer 
contributions. 

The primary difference between the above examples 
appears to be the conductivity at the electrodes. Since the 
equilibrium conductivity has such a strong dependence on 
temperature, the electrode temperature is a critical 
parameter in controlling the sheath losses. A very accurate 
comparison can only be made using a plot generated from the 
correct charge density and on precise knowledge of anode 
temperatures. 
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losses 



It is common practice to predict the 
according to the empirical formula 

V = A + B6 ' (46) 

where d is the boundary layer thickness and A and B are 
constants to be determined for a particular experiment. 
Such a scheme might be useful here if A is the sheath drop 
and B is the boundary layer drop. Using the comparisons 
above and from Appendix A for the boundary layer drop, the 
formula wculd look like 

Sonju and Teno V = 18.8 + 1190(5 

Argyropoulos et al V = 6.1 + 1110(5 
where § is measured in meters and V is measured in volts. 

B. NUMERICAL CONCLUSIONS 

The numerical technique is discussed in detail in 
Section III. Generally it consisted of a finite difference 
scheme for the solution of a system of three (later four) 
coupled, non-linear, partial, second-order differential 
equations. The line iterative method of solution was used. 
This method allowed a smaller storage requirement, and a 
choice of a 51x51 array balanced the computer time to 
converge with a mesh size sufficient to show the detail of a 
two-dimensional sheath. Convergence required four to six 
hours of computer time on the IBM 360/67 computer for each 
potential and/or Kail parameter chosen. The addition of 
Joule heating reduced the convergence time to one hour 
because of the damping effect of the temperature term near 
the electrode. 

The use of successive over-relaxation to reduce the 
convergence time met with only marginal success. Because of 
numerical i nstabilites, an over-relaxation parameter greater 
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than 1.1 usually caused the process to diverge. 
Consequently, a value of 1.05 was used resulting in less 
than a 10% improvement in time. 

In its present form, the compu 
outlined in this work is limited as to 
parameters describing the system. For 
resulted for non-dimensional potentials e 
parameters greater than 2, and Joule 
greater than 10*. It can be observed fro 
Hall parameter weights the relative’ 
derivatives compared to the ccoss-cha 
giving rise to numerical instabilities fo 
values of . The maximum potential attai 
limited by the size of the array since a 
obtained for a 51x51 array and only 29 v;a 
41x41 array. 

C. RECOMHENCATIONS FOR FURTHER STUDY 

Since the principal limitations to this work arises 
because cf the large storage and time required for 
calculations, much more could be accomplished with improved 
numerical techniques. A search by the author for a packaged 
subroutine capable of solving the computational matrix more 
efficiently than the one used here proved futile. Improved 
techniques are needed which will accomodate larger fields of 
computation. 

The validity of these results is limited by the fact 
that the current density decreases in two dimensions rather 
than three. A three-dimensional scheme, admittedly a major 
task to develop, could render more realistic information on 
the effects of the sheath and on other electrode phenomena. 
One possible scheme is to model a three-dimensional system 
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by using an axi-sy mmetric cylindrical geometry for the point 
electrodes. 

Additional physical effects involving changes in 
boundary conditions are possible. For example, boundary 
conditions could change the anode to a cathode b/ changing 
the sign cf the charge and making the electrode an emitter 
of electrons.[ Ref . 27, p. 106] The addition of convection 
into the system of equations (should it be deemed 
appropriate) is expected to accelerate convergence because 
it would ease the requirements of current density decreasing 
away from the electrode. 

Other effects such as neutral and ion heating may be 
incorporated into the scheme to investigate thermal 
instabilities known to exist for such a case. Conduction of 
heat in the solid walls, ion slip, and frozen flow may be 
added. It is anticipated that frozen flow would decrease the 
importance of the boundary layer drop and increase the 
importance of the sheath drop. 

With each additional effect the model rises in 
complexity and the numerical scheme becomes more unwieldy; 
therefore, only the principal effects should be reflected in 
the description of the sheath. 
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APPENDIX A 



A SIWPLIIIPD TECHNIgUE FOR D^TERHING THE BOUNDARY LAYER 
VOLTAGE DROP IN MHD GENERATORS 

As a complement to the two-^dimensional sheath prcblaia, 
and in order to assist in understanding the overall loss 
problem, a one-dimensional boundary layer analysis was 
investigated. The result is a method for determining tke 
voltage drop across the thermal boundary layer which is 
simpler than most methods in existence today. 

The technique balances sophistocation with ease of 
application. Because of its lack of rigor, many simplifying 
assumptions have been made which render the results 
inaccurate for certain cases. These assumptions are 
detailed in the following paragraphs, and experimental 
comparisons are shown. 

1 . Geometr 

Figure A1 is a schematic of the voltage drops across a 
typical channel. the geometry assumes opposing electrodes, 
but does not preclude the possibility of diagonal 
construction. This analysis looks at only the ohmic portion 
of the voltage profile and assumes that the sheath losses 
can and must be determined separately. The boundary layer 
may be fully developed and uses the wall conditions as one 
of the boundaries for ohmic analysis. 

2. Boundary LaiSE Drop 

Reconstructing and modifying a formulation by Ecsa[35], 
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Figure A. 1 Three representations of voltage losses in the MHD generator, 
a) sketch shoving regions of loss inechanisms, b) electrical 
analog, c) voltage profile. 



across the channel can be 



the average resistivity ((7~^) 
written as 



^ 6 



- - 1 . _ 1 , 1 , _ 1 1 

^ D / a a D 



D-6 



/ (^) dy + / 



D 



0 
a 



+ 






D-6 



(f ) dy 



(A,1) 



where 5 is the thermal boundary layer thickness and D is the 

channel width between opposing electrodes. Assuming 

constant conductivity in the free stream region, symmetric 

boundary layers at opposing walls, and non-diraensionalizing 

the integration variable, the equation becomes 

1 

1 + 2 |- / d(|^) - 2 I- (^-2) 

u o 0 



(a ) - - 



c L 



D 



0 



Eosa [Eef. 32, p. 75] assumed a plasma where the gas 
properties vary only across the channel, and averaged Ohm's 
law. The resulting equations for current are 



= a (Ey + U X B) - I 



(A. 3) 



J = a E + 3d,, 
X X y 



(A. 4) 



_ 14/32 

where G= 0{~q~) and all rerras have been averaged 

across the channel (y-direction) . The x-direction is the 
direction of flow. Elirainaring J , assuming Q is a 

X 

constant across the channel, solve for E to get 

y 



UB + -^ (G + 3^^) = 3E 
o ^ 



(A. 5) 



The actual output voltage is 
D 



V 



out 



= / Edy = ED = UBD + (G+3^) + ^Ej,D 



(A. 6) 
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The intercept of the voltage profile with 
(see figure A1) is given by 




the channel wall 



(A. 7) 



Taking the difference and substituting for G the voltage 
loss is 

Vioss = ■ (1 + B )[a^(a - 1] (A. 8) 

c 

This loss is Ohmic and is attributed to the two boundary 
layer temperature drops. For a turbulent boundary layer the 
first term will be small compared to the second 

Vfoi = (1 + 3 ) tcr^(cr ) - 1] (A. 9) 

c 



3 . Co nducti vi ty Ex pression 

a 

An expression will now be derived for to be used in 
Eg. A. 9. Hurwitr, Kilb, and Sutton[36] have shown that for 
cases where the current is entirely cross-channel (Faraday 
mode) , the conductivity can indeed be treated as a scalar 
since for most such channels the magnetic field is 

approximately uniform between the electrodes. For the case 
of scalar conductivity we can write[ 26 ], [ 37 ] 

0 = ^ ^ (A. 10} 

k 

The subscript k denotes the various species present in the 

plasma, and x = n /n . For non-reacting plasmas the sum of 
k k n 

the products x Q will be relatively independent of 

k k 

temperature. We will assume this to be approximately true 
for reacting plasmas. Taking the ratio of the conductivity 
at a point to that at free stream. Eg. A. 10 leads to 
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c 

a 



n c n 
eo e n 



n - 
e eo nO 



(A. 11) 



Assuming a locally Maxwellian velocity distribution, 
guasi-equilibriura (allowing the use of the Saha relation) , 
and constant pressure across the boundary layer. Eg. A. 1 1 
becomes 



_£ = a-V4 a(l/e-l) 

a ^ 



(A. 12) 



Ct 

where CL = and 0 = T/T . 

2kT„ o 



4. Constant Spe ci fic Heat 



At this point a description of the turbulent boundary 
layer is needed. For the purpose of this analysis we will 
use the 1/7th-power law and Reynold's anaiogy[ 38 ] to 
describe the temperature field. An extensive and detailed 
study of the heat transfer problem related to these high 
temperatures was conducted by Brira[39], but in the interest 
of simplicity the Prandtl Number is assumed to be one. For 



constant C , 
P 



Z = 
6 



T - T 



w. 



‘T - T ' 
o w 



0 - 0 , 



’V7. 7 



1 - 



6w 



(A. 13) 



so 



' - <fw> 



d0, 



(A. 14) 



Now combining Egs. A. 2 and A. 12 and changing variables 
according to Eg. A. 14, the result is 



(o = 



o 



e 



w 



w 



(A. 15) 



Putting Eg. A. 15 into A. 9 for the average resistivity gives 
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(A. 16} 



^bl = 3^)6/a^ 



7_^ / 0-^/4 de 



(1 - e ) 0 

^ w w 



w 



- 1 



5 • S pec ifi c Heat 



is 



For a plasma with a non-constant C the 1/7th-power law 

P 




h - h 




7 



(A. 17) 



A change of variables would then give 

L»'o-V J 

Obviously an enthalpy table or equation is needed to solve 
Eg. A. 18. However, Mollier diagrams are obtainable, and 
once the enthalpy relation is known, the method shown in 
this analysis retains its simplicity. 



6 • E Fu ncti ons an d Voltage Drop Para met ers 

Figure A2 is a graph of the conductivity from Eg- A. 11 
for Toluene and pure oxygen with 2% cesium seed at a core 
temperature of 2600OK. Plotted also is the con ductii'ity 
from a more soph is tocated program provided by AVCO Everett 
which also assumes guasi-eguilibrium, but does not make the 
other simplifications done above. The agreement at 260QOK 
results from the "forced" condition in the simplified 
program . 

An examination of Egs. A. 16 and A. 18 shows that the 
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Figure A. 2 Conductivity for stochiometr ic nixture of Toluene and Oxygen with two percent 
Cesium seed at 0. 3 atmospheres . Base temperature = 2600^K. 



effect cf changing the integration variable was to average 
the resistivity across the boundary layer in "teiaperatarc 
space" instead of "geometry space". In each case this change 
generated a weighting function which weights the 
contribution of the resistivity to the integral according to 
the temperature. Other weighting functions derived here 
are : 



for constant C : 
P 






for variable C : 

P 

Other descriptions 
Useful descriptions 
available in Cramer 



W = 7C T 

h p o(h„ -h^) 7 

are possible, and often desirable, 
for laminar boundary layers are 
and Pai[40] and Kerrebrock[ 4 1 ] . 



Figure A3 shows the effect that the weighting function 
has on resistivity. The resistivity (reciprocal of 
conductivity) is plotted first from the AVCO program and the 
approximate value from Eg. A. 12. Next, the reacting plasraa 
weighting function times the differential Ao is plotted. 
Lastly, the product of the first two graphs shov/s that the 
effect is to reduce the error considerably, by weighting 
more heavily those values near the core. In fact the 
average resistivity from Fig. A3 (area under the right hand 
graph) is 0.0254 for the AVCO program, and 0.0244 for the 
simplified program. Obviously, the effect of the boundary 
layer shape is to smooth out the errors of the approximate 
resistivity calculations. 



From Egs. A. 16 and A. 18 a non-dimensional boundary layer 
voltage parameter can be identified. 



4 - 



bl 






(A. 19) 



so that 
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Figure A. 3 Effect of weighting function on two programs for resistivity of Toluene/Oxygen/Cesium 
plasma using l/7th power law and varying C . 



\ 



W{0) d6 - 1 



(A. 20) 



4-k. = O'’/'’ 

bl ' 

0 

w 

where W (0) is the weighting function appropriate to the 

boundary layer description. Figures A4a and A4b are plots 

of 4> vs. wall temperature and core temperature for a 
bl 

constant C gas and Cesium seed. The former uses the 
P 

1/7th-power law while the latter uses a 1/8th-power law. As 
might be expected the 1/8th-power law predicts less of a 
voltage loss because the region of lower temperature plasaia 
is thinner. Figures A5a and A5b are similar graphs for the 
combustion products of Toluene in pure oxygen and 2/’? cesium 
seed added. 

7. E xpe rimen ta l Comparisons 

As with all experimental predictions the accuracy of the 
results is dependent upon the accuracy of the input data as 
well as the ability of the model to represent the actual 
situation. The following exaraples illustrate the 
convenience of this method, but at the same time show the 
importance of an adequate model. 

Sonju and Teno[34J report the results of an experimental 
run on the Viking I generator. Fig. A6 is a reproduction of 
the voltage profile for one run. Because of burner 
innef f iciencies they found the conductivity of the 
Toluene/oxy gen/cesiura plasma to be between 12 and 16 mhos/m, 
about half the 33 mhos/m that would have been predicted for 
complete combust ion[ 42 ]. Using Q =1.8, J = 1.6x10'^ 

y ■ 

amps/m2, = 12-16 mhos/m, S~ *1® {half the channel widrh 

for fully developed flow), T =2500<>K, T =180QO, The 

core wall 

voltage drop parameter is 0.282. This equates to a voltage 
loss per boundary layer of 119 volts forCT=l6, and 159 volts 
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NON-DIMENSlO.tMAL POTENTIAL DROP <P 



G. 




non-dimensional potential drop <(> 

b. 



Figure A. 4 Non-dimensional voltage drop (}> versus wall temperature and 
free stream temperature for Cesium seed, constant C^. 
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NON- DIMENSIONAL POTENTIAL DROP <|> 



a . 




b. 



Figvire A. 5 Non-dimensional voltage drop <$> versus wall temperatxare and 
free stream temperature for Cesium seed, varying C^. 
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TRANSVERSE VOLTAGE- VOLTS 



300 




Figure A-.6 Observed transverse voltage distribution. 
(After Sonju and Teno ) 
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for (7^-12. Extending the voltage profile to the channel wall 

in Fig. A6 shows a V of about 150 volts. 

bl 

Another example comes from Argyropoulos et al[33] who 
use a sophistoca ted computer technique to compare theory 
with experiments. They arrived at a boundary layer 
potential drop at the anode of 101 volts for an exj^riment 
on the AVCO-APL channel. Using the same data, namely 

T = 2600OK, T = 2000OK, Q = 1.02, (5 =60 mm, a 

core wall ^ 

chemistry of toluene, oxygen and cesium, and the 1/7th-power 

law. Fig. A. 5 gives a non-dimensional potential drop of 
0.24. This corresponds to a predicted drop of 110.5 volts. 

The next example shows the importance of selecting a 

proper model for the boundary layer. Fig. A7 is from data 

of Kessler and Eustis[6J. The parameters for the experiment 

were: I = 2.75 amps, electrode surface area = 3.68 cm2, 

T = 16850K, T =27000K, (^ = 10 mhos/m, /3=1.5, The gas 
elec core ^ 

was Ethanol in oxygen with ^% KOH added to increase 

conductivity. Nitrogen was introduced for cooling such that 

the N /O ratio was 0.5. Using the 1/7th-power law 4> was 
2 2 bl 

found to be 0.290. Kessler predicted the boundary layer 

thickness to be 1.2 cm. The calculated voltage drop is then 

8.5 volts falling well below the 46 volts shown in Fig. A7. 

The difference is too great to be attributed to sheath 

phenomena alone. 

The geometry of this last case is the cause of the 
discrepency. The channel depth is 3.1 cm while the width 
between opposing electrodes is 10 cm. This means that the 
flow does not behave as between two flat plates since the 
sidewalls have a strong effect on the electrode wall 
boundary layer. The model used in this comparison is 
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POTENTIAL RELATIVE TO WALL CENTER 

VOLTS 



Figure A. 7 Sidewall probe potential, <^re temperature = 2700^K. 
(After Kessler and Eustis ) 
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obviously inappropriate. Schlichting (Ref. 38, p. 575) 
contains information on the problem of wall interference in 
non-circular cross-sections and points out the increase in 
friction and boundary layer size. A gross, but closer 
approximation for this procedure is to assume fully 
developed flow (as was suggested by Reseck[43] for this same 
channel). This results in a voltage drop prediction of 35.5 
volts. Kessler attributes 10 volts to the sheath drop for 
the cathode. When combined with the 35.5 volts predicted 
above, the entire drop is accounted for more fully. 

8. Conclusions 

Equation A. 20 represents a useful means of determining 
thermal boundary layer voltage losses as long as the 
following conditions are met; 

1) A function of enthalpy vs. temperature is available, 

or, for the case of Eg. A. 16, the gas has a fairly 

constant C in the range of temoeratures of the channel- 
F 

2) the core conductivity is known or can be estimated with 
some degree of accuracy, 

3) Hall parameter, core temperature and wall temperature 
are known, and 

4) a suitable boundary layer description is available. 

The simplicity of this approach makes it desirable for 
integration with a computer program which predicts the 
performance of an HHD generator. Such programs usually 
already contain free stream and wall temperatures, core 
conductivity, and Hall parameters. Since these equations 
are one-dimensional they allow calculations at stations 
along the channel, eliminating the need to determine 
properties at points in a two-dimensional field. A sample 
program applying this method at one station is given in the 
section "Computer Programs." 
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APPENDIX 3 



PROOF OF TiiE N ON -EXISTENCE OF A ONE-DIHENSIOEA L SH^Til 



A first attempt at solving a complicated set of 
equations is usually to reduce the problem to uiie dimension. 
For the equations which describe the sheath this procedure 
seems reasonable due to the small thickness of the sheath 
relative to the electrode dimensions. Such an attempt for a 
constant temperature plasma with no ionization or 
recombination will fail because of the nonexistence of such 
a solution. P.esults published in Ref. 21 are summarized 
here along with significant conclusions which were 
instrumental in shaping the course of the tv.o-di mensionai 
work. 

In order to simplify the analysis a change of variables 
is made. 



The following manipulations are made on Sgs. 26 and 27 in 
the text: 

1) Magnetic field effects are dropped, 

2) The equations are considered in only one-dimension, 

3) The equations are added to each other and subtracted, 

4) The change of variables in Eq. B. 1 are applied. 

The results, including Eq. 8 are 



SOLUTION 



r 




E = n . + 

^ 1 e 



(B.1) 




0 



(B.3) 



(B.2) 
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0 



(B.4) 



dy 



(C 



dy^ 



d^r 



dy' 



If Egs. B.3 and B.4 above are integrated and E = - d4>/dy 
is introduced then 



^ = r 

dy 



(B.5) 



- TE + ^ = C, 

dy 1 



(B.6) 



- CE + 



dy 



= C, 



(B.7) 



It is not difficult to show that, in the absence of a net 
ion current 

Cl = -C2 = Jy/(eD^) 

where J is the conventional current flux and D is the 

y e 

electron diffusion coefficient. Obviously, since J exists, 

y 

C and C are non-zero constants. 

1 2 



Because the ambipolar region is guasi-neutral, F 
approaches zero away from the anode, reaching a value of 
zero at the plasma proper. Thus, as y increases: 

r approaches 0 
^ approaches ^ 

E approaches 

Eventually a distance y = y would be reached so that for 

o 



every y > y 



jrEj < C^/2 



(B.9) 



Hence 
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1 



(B. 10) 



dy 



> 



Cj^/2 



for all y > Yq 



5(y) i C{y^) + (Cj^/2) (y -Yq) 



(B. 11) 



which does not approach a constant as required. 
Furthermore, since 



dy 



-C^ + CE 



(B. 12) 



^ > -Cl + + Ci/2 X (y-y^)] E , (S-13) 



for y > y , which eventually grows at least 
c 

'^1 2 

E ^ (y - Yo' 

and does not approach zero as required. 



as fast as 

(B. 14) 



These obvious inconsistencies arise because of the 

constant C . The necessary requirement for a solution is 

1 

that the current density decrease away from the anode giving 
rise to current constrictions at the surface. Under these 

conditions, C would decrease at an appropriate rate so that 

1 

r approaches zero and ^ approaches ^ . Mechanisms which 

allow a proper variation of C are diffusion in two or three 

1 

dimensions and ionization/recorabination. The former negates 
the applicability of a one-dimensional Cartesian solution 
while the latter eliminates the frozen flow assumption. 
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APPENDIX C 



SIMI'NS-ION AL AND FRACTIONAL AN A LYSES OF THE SHEATH EQOA TI ONS 

Eguations 1—6 in the text are cumbersome in tliat they 
contain many physical parameters which in their present form 
must be carried along as the equations are operated on. It 
is useful to conduct a dimensional analysis or these 
equations in an attempt to collect the coeff icieiits and 
variables in such a way as to simplify the appearance of the 
equations and better understand their nature. Also, before 
introducing the problem into the computer it is imperative 
to know the sizes of the quantities that are being dealt 
with. For this latter reason, a fractional analysis showing 
the relative significance of certain physical effects is 
also given in this appendix. 

The method outlined by Langhaar[44] is utilized here to 

determine the non-dimensional parameters relevant to the 

sheath equations. The significant parameters whicn appear 

in Egs. 1 and 3 are J, <|) , n , y, e , f , D , ^ , R, and B, 

s o s s ^ 

where y represents any coordinate length. An expansion and 

manipulation of Eg. 1 with Egs. 3 and 5 will show that the 

terms D and U can be replaced by kT . Since k and T 
S s o 

always appear together, they can be expressed as a single 

tern. Additionally, the Hall parameter p is already 

non-dimensional and will not be included in this analysis. 

This leaves n , y, e, £ , and kT . In the 3KS system 

s o o 

these parameters are combinations of the units Oj. -mas^ (H) , 

length (L), charge (g) and time (t) . The matrix showing the 
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coefficients of these units for each parameter is: 





4> 


n 


L 


e 


6 

0 


_ 0^ 


M 


1 


0 


0 


0 


- 1 


1 


I 


2 


-3 


1 


0 


-3 


2 


q 


- 1 


0 


0 


1 


2 


0 


1 

t 1 


L:iJ 


0 


0 


0 


2 


-2 



At 



irst 


glance 


it might 


be expected 


that the 


= 2 


dimensionless parameters (nurabe 


r of par 


the 


number 


of units 


necessary to 


describe 


er. 


closer 


examination 


shows that 


the t row 


rly i 


ndependent of the W 


row. In fact 


the rank 



matrix is three ijot four. Thus there are 6-3 = 3 
independent dimensionless parameters (number of parameters 
minus the rank of the matrix) which describe the system of 
equations. This result would also have arisen if k and T 

o 

had been considered separately. 



Disregarding the t row, the homogeneous linear algebraic 
equations whose coefficients are the numbers in the rows of 
the dimensicral matrix are 



2a^ - 3a2 ^3 ~ ^^5 ^ 

“^1 ^4 ® 

Solving for a , a , and a in terms of the other three gives 
4 5 6 
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ai + 6a2 



^4 

^5 



- 2a, 



-3a2 + aj 
-«1 - 3aj + aj 



(C.2) 



Now setting the values a =1, and a =a =0, the result is; 

1 2 3 

a =1, a =0, and a =- 1 . 

4 5 6 



Continuing by assigning a value of one successively to 

a and a with the other two of that set equal to zero, the 
2 3 

solution matrix is generated. 



({) n L e 6 kTp 

o 









i 








A 


1 


0 i 




0 


-1 


A 

n 


0 i 

j 


1 

i- J 


0 1 6 
1 


-3 


-3 


A 

y 


0 1 


L.!.j 


1 i -2 

i J 


1 


1 ' 



The resulting non-dimensional parameters are then 







^ f . 3 

" = 'krr:’ " 

o o 



e kT 
o o 



(C. 3) 



These terms are all independent since n and y appear in 
only one of each. In numerical form they look like 

$ = (j>/({)^ = 1.16x10^ (}) 

- / Q in-12 m "3 „ (C.4) 

n = n/n = 9.26x10 T n 
o o 

J = y/L = 4.76x10^ y 
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2 . 



E st imate of the Sheath Size 



The sheath size can be estimated by a techni< 3 ue called 
fractional analysis or approximation theory[45]. Eq. 3 is 

„2 



•(J) = -(n^ 



(C.5) 



Now "characteristic" parameters <j) , ^ , and n are factored 

os o 

out of the equation for each variable. 



‘J’o '- 2 '' 

TT ^ ^ 



n e ^ 
o i'' . 

— — (n. - n^) 
e 1 e 
o 



(C.6) 



where <}) and n are the maximum value of the potential and 
o o . 

charge concentration respectively, and X is a sheath 

s 

cha ractersitic length. Since in the sheath all the 
variables are of order unity, for the eg-uation to be 

n„eA^ 

consistent the term ~r must also be of order unity. 



Setting this term equal to one and solving for A , ti 



le 



result is 



X -Q Q (C.7) 
s n^e 

This gives an order of magnitude for the characteristic 
sheath length. For example, for a sheath drop of 10 volts 
in a plasma ionized to electrons/tu^ the characteristic 
sheath length is 2.4xl0~5m. 



3 . 



Characte ris t 
From Eg. 10 



ic Current from Conduction 



the 

J 



e 



conduction current 



is 



(C.S) 
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Working 

current 



in one-dimension and identifying a characteristic 

J , a fractional analysis gives the equation 
o 



J A 
o s 

(b u en 

o 



(C.9) 



Notice that in order to be consistent the char a ct eristic 

length in the potential derivative roust again be A , the 

s 

sheath length. Solving for J when n = jL/ = 1 

o o e 

m2 (volt-sec) -1 , and (b = 10, gives J = 6.8x10'» amps/ci2. xhe 

o o 

value for jj comes from Eg. 16 for electrons colliding with 
e 

CO molecules at IOOQok. 

2 

4. Diffusion Length 



The potential is a monotonic positive function which is 
suitable for approximation theory, however the charge 
density profiles have many inflection points in the sheath 
and are not well suited. Nevertheless, it is possible to 
consider profiles beyond the sheath where characteristic 
lengths can be defined. 



If diffusion 
production then 



is to play a significant 



J = y 
o e 



8n 

8y 



role in current 
(C. 10) 



Again, using fractional analysis, the following 
approximation is made: 



Solving for 

2.03x10-7. 

played by th 
within the s 



o a 



n^y^kT^ 
o e e 



1 with the same parameters as 



d 

Since L <A it can be concl 
d s 

€ mechanism of diffusion s 
heath. 



(C. 1 1) 

above gives L = 
d 

uded that any role 
hould be visible 
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5- lonizat ion /Rec om bina tiog L en gth 



From Hinnov and Hirschberg[ 46 ] the three-body 
recombination rate of electrons and ions is given by 

dn. 



e 2 

■ 3 T- ■ = an 

dt e 



(C.12) 



where 



c c 1 n~33 /kT. -4.5 
a = 5.6x10 ' ( — ) n 

e e 






(C.13) 



Using the technigue outlined previously, 

_ c: ,kT>-4.5 3-3 

— jir- = 5.6x10 ( — ) n n 

dt e o e 



(C. 14) 



where T is the characteristic recombination time. Let the 
r 

term 5 . 6x 1 0~39 -4 «. Sfiz Le of the order of unitv and 



solve for T . 

r 



/c r in”33. ,kT.-4.5, 2v~l 

= (5.6x10 ) (™) (n^ ), 



(C.15) 



At T=1000®K and n = 10^®m~3^ the characteristic tirae 

o 

becomes 2.9x10-2sec. 



In order 
characteristic 
since this is 
recombination. 



to change the characteristic time to a 
length the mean thermal velocity is used 
the transport mechanism for ionization and 
Thus 

/8kT ' 

(C.16) 



L = T c = T 
r r r 



8kT' 

e 



irm 



Using the same conditions as above L = 5.7x102 m. Since 

r 

L >>X it can be concluded that ionization and recombination 
r s 
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do not have time to occur within the sheath for the 
conditions used in this analysis. Since the results of this 
work show that the ambipolar region is of the same order as 
the sheath length, ionization/recombination is not expected 
to play a part in that region either. 

6. Convection Time 



At atmospheric pressure and temperatures around lOOD^K 
sonic velocity is of the order of 10^ m/sec. This is the 
order of velocity of the flows of MHD generators. A rough 
characteristic time relevant to the sheath could be 
established by comparing the flovf velocity with the size of 
the sheath. 



T 



c 



X /U 
s c 



(C, 17) 



Using (j) =10 volts and n 

o o 



1018 tn-3, T = 2.4 x10-8sgc. 
c 



A characteristic sheath time can be defined from 



= X /c 
s s' 



(C. 18) 



for comparison with T • Using representative values for X. 

c s 

and c, T = 1.2x10-io sec. Thus, any effects occuring from 
s 

convective phenomena such as mixing and boundary layer 
effects do not have sufficient time to affect the sheath. 
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COMPUTER PROGRAMS 



The principal program used in this v/ork is the one that 
calculates the potential and charge density distributions 
for the sheath and ambipolar regions. This program is 
explained in Part 1 of this section and reproduced there. 

Auxiliary programs are those that calculate the current 
distribution and produce the outputs used in this work. 
Part 2 shows the auxiliary programs. 

Part 3 is a reproduction of the program used in the 
calculation of the boundary layer in Appendix A. A sample 
enthalpy table is provided in data format for the 
Toluene/oxy gen/cesium mixture. 

^ ♦ Sheath Pro gram 

The sheath program begins v;ith an assumed solution and 
operates on the field of parameters to replace the assumed 
solution with a more accurate one. Theoretically, the more 
iterations that are used, the more accurate the result is. 
An example of an initial assumed solution is P{T,J)=0, 
G(I,J)=BZ, SI (I, J)=BZ, and TH(I,J)=1.0 throughout the field, 
(see the program comment cards for definitions) These three 
arrays are operated on and a more accurate solution is 
stored for later analysis or further iterations. 

The following pages are a reprint of this program, and 
an explanation of the parameters and subroutines is 
included . 
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OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 



SHEATH COMPUTATION PROGRAM 






THIS PROGRAM IS DESIGNED TO CALCULATE A TWO-DIMENSIONAL 
SHEATH FIELD GIVEN THE FOLLOWING INITIAL PARAMETERS: 

P POTENTIAL DISTRIBUTION 

G ELECTRON DENSITY DISTRIBUTION 

SI ION DENSITY DISTRIBUTION 
TH ELECTRON TEMPERATURE DISTRIBUTION 
BET HALL PARAMETER 
GAM JOULE HEATING PARAMETER 
ZL STEP SIZE multiple 
PHI ELECTRODE POTENTIAL 

W SUCCESSIVE OVER-RELAXATION PARAMETER 

BZ CHARGE DENSITY BOUNDARY VALUE 
MAX NUMBER OF ITERATIONS PER CYCLE 
MM NUMBER OF CYCLES 

JMIN MINIMUM PRINTOUT PARAMETER (BEGINS PRINTOUT AT 
END OF ITERATION JMIN WITHIN EACH CYCLE) 

NRCWS ARRAY WIDTH 
NCGLS ARRAY LENGTH 
NCEM ELECTRODE LOCATION ON WALL 
THIS ROUTINE USES ASSUMED VALUES F CR ARRAYS P,G,SI AND TH 
AND UPDATES THESE VALUES. THE PROGRAM USES THE FOLLOWING 
SUBROUTINES: 

MAIN. . .CONTROLS INPUT AND OUTPUT AND DEFINES MOST 
CONSTANTS AND PARAMETERS 

LOAD.. .LOADS THE COMPUTATIONAL MATRIX AND CONTRCLS 
THE SWEEP MODE 

OERV. . .DERIVES INTERMEDIATE DERIVATIVES FOR USE IN 
LOAD 

RCTAT. .ROTATES ARRAY FOR CROS S- COG BO I NATE SWEEPS 
TMP. .. .SOLVES ENERGY EQUATION FOR USE IN LOAD 
WRITER.. PR INTER OUTPUT SUBROUTINE 

OGELB..I3M PACKAGED SUBROUTINE FOR THE SOLUTION OF 

NON-SYMMETkIC SIMULTANEOUS ALGEBRAIC EQUATIONS WITH 
BANDED STRUCTURE. 

IMPLICIT REAL«8( A-H,0-Z) 

COMMGN/BLl/ P( 51 ,51)tG(5L,5i)tSI(51,51) 
CCMMON/BL2/3ET,GAM,H ,PHI » BZ t NCOLS , NRCWS, NCEN , MAX 
COMMON/ BL 3 / BM AX , W , JM I N , J MAX 

CCMM0N76L4/PP( 51 ) ,PPP (51 ) , GXP( 51 ) tGXPP(51)»SP(51),SPP{ 
51 ) 

COMMON/BL5/NLFT,NRIT 
C0MMCN73L6/TH(51 ,51 ) , DTX( 51 ) 

C READ CONSTANTS 

READ( 5,200) PH I , ZL , B E T , W ,GAM , BZ 
READ( 5,3DC) MAX, MM, JMIN, NCEN 

WRITE(6,201) PHI.ZL, 3 ET , W , GAM , B Z , MAX , MM , J M I N , NC EN 

JMAX=MAX 

BMAX=1 .00-08 

TEMP=2.00 03 

NPCHS=51 

NCCLS=51 

H=ZL/OFLOAT(NROWS-l ) 

PHI=PHI*1 .1 609D 04/TEMP 
00 1 J=l, NCOLS 
OTX( J )=0.0 
PP( J)=0.0 
6XF{ J)=0.0 
1 SP(J)=0.0 

DO 11 J=l, NCOLS 
DO 11 I=1,NR0WS 
TH( I , J) =1 .0 00 
11 CONTINUE 
C READ STORED ARRAYS 
READ(2) P 
R£A0(2) G 
READ(2) SI 
REWIND 2 

C OPERATE CN ARRAYS 
DO 2 KK-l,MM 
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CALL LCAD 
2 CONTINUE 

C OUTPUT AND STORE RESULTS 
WRITE(2) P 
WRITE(2) G 
WRITE(2) SI 
WRITE! 6, 100 ) 

100 FORMAT (3Xt • DATA HAS BEEN READ INTO FILE*) 

200 FCRMAT(6D10.0) 

300 F0PMAT(4I4) 

201 F0FMAT(3X,6012.3,4I4) 

STCP 

END 



C 

C 



C 



C 

C 



SUEROLTINE LOAD 
IMPLICIT REAL*E( A-H 7 O-Z) 

REAL*4 EPS 

C0MM0N/BL1/P( 51, 51 ) ,G (51, 51), SI (51, 51) 

C0MM0N/BL2/BET ,GAM,H ,PHI , BZ , NCOLS , NROWS , NCE N , MA X 
C CMMGN/ BL 3 / B.M AX , W , J M I N , JM AX 

C0MH0N/3L4/ FP( 51 ) ,PPP(51),GXP(51),GXPP(51),SP(51),SPP( 
51) 

CCMM0N/BL5/NLFT,NRIT 
COMMON/BLS/TH( 51 ,51) ,DTX(51) 

DIMENSION T(153),A(10,153),B(1505) 

EQUIVALENCE ( A ( 1 , 1 ) , B ( 1 ) ) 

HS=1 .CD 00/ (H*H) 

HR=.5C OC/H 
H4=HS/4.D 00 
BETR=BET-HR/Z.D 00 
NCCM=NC0LS-1 
NFLAG=1 
NCCLS3=10 
NRC«S3=NR0WS*3 
MU0=4 



MLD=5 

BEGIN ITERATIONS 
DO 99 JCNT=1,MAX 
ND = 2 

ROTATION parameters 

I N = NF LAG- (NFL AG/4) *4 
IF( IN.EQ.0.CR.IN.EQ.3) ND=1 
IF(GAM.LT. 1 .D 00) GO TO 30 
CALL TMP(ND) 

SWEEP IN I DIRECTION 
30 DC 1 IX--1,NR0WS 
I = IX 

IF( IN.EQ.0.0R.IN.EQ.2) I=NR0HS-IX+1 
DC 90 M=1,NR0WS3 
DO 90 L=1,NC0LS3 
A{ L, M)=0.0 
90 CONTINUE 

IF(I .NE.l .AND.I.NE.NROWS) CALL DERV(I,2) 
DO 2 J=l, NCOLS 

INTRODUCE TEMPERATURE PARAMETERS 
TT=TH( I , J ) 

THS=TT*TT*HS 

LOAD CCMFUTATI CNAL MATRIX 
DO 2 K=l,3 
M = 3*( J-1) +K 



GO TO (50,51 ) ,ND 

50 IF( J .EQ.l . AND. I. EQ.NCEN) GOTO (76, 21, 21), K 
IF(J. EC. NCOLS ) GO TO (21, 19, 19), K 

IF( I .FQ.l ) GO TO 91 
IF( I .EQ.NROWS ) GO TO (96, 97, 98), K 
IF(J.EQ.l) GO TO (40, 19, 19), K 
DTY = TH(I,Ji-1)-TH(I,J-1) 

TS = .5D 00*TT*0TY=i‘H4-BET’i'-TT«HR=i'( PP( J) +.5D OO-'i'DTX ( J ) ) 
GO TO { 5, 6, 7) , K 

51 IF( I .EQ.l. and. J. EQ.NCEN) GO TO (76, 21, 21), K 
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IF ( I . EQ.NROWS ) GO TO 1 
IF(J.EG.l ) GO TO 91 
IF ( J. FQ..MCOLS ) GO TO 86 
IF(I.EQ.i) GO TO (41, 19, 19), K 
DTY=Th( I , J + 1 )-TH( I , J-1 ) 

TS = .5C 00=;'TT*DTY*H4 
GO TO (5,6, 7) , K 

C INTERNAL FIELD COEFFICIENTS 

5 A( 7,M)=-1.D 00 
A(8,M)=1.D 00 
A(3,M)=HS 
A(9,M)=HS 
A(6,M)=-4.0 00=^HS 
T(N)=-PPP( J ) 

GO TO 2 

6 A( 3, H)=THS-TS 
A(6,M)=-4.D C0«THS 
A(9, M ) = THS + TS 

A( 2,M ) = TT-'M G( I , J f I ) - G (I , J- 1) )’'-H4-.5D 00*G( I , J ) * DT Y^!'H4- 
BET^HR-'.'^ ( TT=t=GXP( J )-.5D 00=:'G( I , J )*DTX( J ) ) 

A( 8, M)=-A( 2 ,M) 

T(P) = TT«(-TT*GXPP( J) <-GXP( J )^"PP( J)- 

.5D 00*DTX( J l*GXP( J)+G( I , J ) =;»( G( I , J ) -S I ( I , J ) ) ) - 
.5D 00^=G( I, J)-OTX( J)=-PP( J)- 
.5D OO^BET^i'OTY^HR^IGI I , J ) '-^P P ( J H- 
TT«GXP( J) ) 

GO TO 2 

7 A(3,M)=HS 
A(6,M)=-4.D OO^HS 
A( 9, M) =HS 

A( l,M)=-( SI ( I , J + D-S I ( I , J-1) )^-H4 
A(7,M)=-A(1 ,M) 

T(P)=SI(I,J)=!MSI(I,J)-G(I,J))-PP(J)-SF(J)-SPP(J) 

GO TO 2 

C BOUNDARY CONDITIONS 
76 T(M)=PHI 
GO TO 85 

91 11=1 
J1 = J 

IF(ND.EQ.l) I1=NR0WS 
IF(ND.EQ.2) J1=MC0LS 
GO TO (92,93,94) ,K 

92 T(N)=P( II, Jl) 

GO TO 85 

93 T(M)=G( II, Jl) 

GC TO 85 

94 T(M) =SI ( n , Jl) 

GO TO 85 

96 T(M)=P(2, J)-P(l, J)+P(NCOM, J) 

GO TO 85 

97 T(P) =G(2, J)-G( 1, J)+G(NCON, J) 

GO TO 85 

98 T(M) = SI(2,J )-SI( 1,JH-SI(NC0M, J) 

GC TO 85 

86 A(3,iM)=-l.D 00 

GO TO (87, 88, 89), K 

87 T(M)=P(I ,2)-P( 1,1) 

GO TO 85 

88 T(M)=G( I ,2)-G( 1,1) 

GO TO 85 

89 T(M)=SI(I,2)-SI( 1,1) 

GO TO 85 

21 T( M) =0.0 
GC TO 85 
19 T(M)=eZ 
GO TO 85 

40 A(6,M)=3.0 OO^BZ 
A(9,M)=-4.D 00*BZ 
A<7,M)=-3.0 00 

T ( M? = -BZ*P? I?3 ) +G( 1,3) +BET«BZ*( P(I + 1,1)-P( 1-1,1)) 

GO TO 2 
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41 T(M)=P( I, J) 

85 A(6, M)=1.D 00 
2 CONTINUE 

C MODIFY MATRIX FOR USE IN SUBROUTINE DGELB 
KMAX=NRCWS3-4 
MB = 0 
MK=6 

DO 200 KA=1,KMAX 
MK = KiK-l 

IF(KA.LT.6) MB=MB+MK 
DC 200 KB=l,10 

IF(KA,LT.6.AND.KB.L6.HK) GO TO 200 

M=10=^(KA-1) +KB 

MA=M-R8 

B( MA) = A(KB, KA) 

200 CONTINUE 
KMM=KMAX+ I 
MK = 0 

DO 202 KA=KFM,NR0WS3 

MB=f'iK<-MB 

MK=MK+1 

MP = 1 0-,VK 

DO 202 K3=l,10 

IFIKB.GT.MPI go to 202 

M=10*(KA-i) fKB 

MA=M-FB 

B(FA)=A(K3,KA) 

202 CONTINUE 
C SOLVE MATRIX 

CALL CGELBI T,B,NR0WS3, 1 , MUD , KLO , EPS , lER) 
IP=I +1 

IFdN.EQ.O.OR.IN.EQ.Z) IP = I-l 
IF(IP.GT.l.ANO.IP.LT.NROWS) CALL DERV(IP,H 
C SUCCESSIVE OVER-RELAXATION FEATURE 
DO 13 J=1,NC0LS 
M = 3^d J~l) 

HINTl = T(Mi-l )-P( I , J) 

HINT2=T(M+2)-G( I , J J 
HINT3 = T(M4-3 J-SK I, J) 

P( I, J ) = P( I , JJ <-N*HINT 1 
G(] , J)=G( I , J ) +W-HINT2 
13 SKI T J)=SI( dJH-W«HlNT3 
I CONTINUE 
C PRINTER OUTPUT 

WRITE(6,103) JCNT 

IF( JCNT.LT. JMiN.OR. JCNT.GT. UMAX) GO TO 43 

WRITE (6, 100) JCNT 

CALL WRITER! Pt NROWS, MCCLS ) 

WRITE(6,101) JCNT 

CALL WRITER (G, NROWS ,NCOLS ) 

WRITE (6, 102) JCNT 

CALL WR I TER ( S I » NROWS f NCOLS ) 

43 IF( IN.EQd.OR. IN.EQ. 3) GO TO 99 
CALL PCTAT(P, NROWS, NCGL5) 

CALL RCTAT(G, NROWS, NCOLS) 

CALL RCTATISI , NROWS, NCOLS) 

CALL RCTAT( TH, NROWS, NCOLS) 

99 NFLAG=NFLAG+i 
RETURN 

100 FORMAT!//, 3X, 'POTENTIAL PLOT', 14,/) 

101 FORMAT !//,3X, • ELECTRON PLOT', 14,/) 

102 F0FMAT!//,3X, • ION PLOT', 14,/) 

103 FORMAT !3X, 14) , , . 

104 FORMAT!// ,3X, • TEMPERATURE PROF I LE* , 14 ,/ ) 

END 
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SUBROUTINE DERV( I ,NF I 
IMPLICIT REAL=s=8( A-H, 0-Z) 

C0MMCN/BL1/P(51,51) ,G(5l,51),SH51,51) 

COMMON /BL2/ BET, GAM, H , PH I , 3 Z , NCOLS , N ROWS , NCEN , M AX 
C0MM0N/BL4/PP (51J,PPP(51) ,GXP( 5U ,GXPP( 51 I , SP( 51) ,SPP( 
51 ) 

C0MMGK/BL6/TH( 51 , 51 ) ,OTX( 51 ) 

LM=NROWS-i 

11 = 1 +) 

12 = 1-1 

IF(I.EO.l) I2=LM 
IF( I .EQ.NROWS ) 11=2 
IF(NF.EC.2) GO TO 4 
C FIRST DERIVATIVES 
H2=.5D 00/H 
DC 1 J=l, NCOLS 

DTX( J) = (TH< II, J) -TH( 12, J) ) =!=H2 
PP{ J) =(P( I 1 , J)-P( 12, J ) )-’--H2 
GXP( J) = (G< I 1,J)-G( 12, J) )^-<H2 

1 SP(J)=(SI( II, J)-SI( 12, J) )*H2 
RETURN 

C SECOND DERIVATIVES 
4 HS = 1 .CD 00/ (H'-'^H) 

DC 2 J=1,NCCLS 
PPP(J)=(P(I1,J)+P(I2,J))^HS 
GXPP( J)=(G( II, J)+G( 12, J) l-'^HS 

2 SPP( Ji=(SI( I1,J)+SU 12, J) )«HS 
RETURN 

END 



SUBROUTINE WRITER! Z , NCOLS,NROWS) 
IMPLICIT REALMS! A-H,0-Z) 
DIMENSION Z (NROWS, NCOLS ) 

Jl = l 
J2 = l 1 

10 IF (J2 .GE. NCOLS) J2 = NC0LS 
DO 15 1=1, NROWS 

15 WRITE(6,220) I , ( Z ( I , J ) , J = J 1 , J2 ) 
WRITE(6,240 ) 

IF(J2. EC. NCOLS) RETURN 
J1=J l+ll 
J2=J2+11 
GO TO 10 

220 FORMAT! • • , I 3 , 1 1 E 1 1 . 4 ) 

240 FORMAT!/) 

END 
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SUBROUTINE ROT AT ( A , N ROWS , NCOLS ) 
IMPLICIT RE^L*8 (A-H,0-Z) 
DIMENSION A (NROWStNCOLS) ,6(51,51) 
DO 1 I=1,NRGWS 
DO 1 J=l, NCOLS 

1 B( I , J)=A( I , J) 

DO 2 i=l,NRCWS 
DC 2 J=l, NCOLS 

2 A( I , J)=B( J, I ) 

RETURN 

END 



SUBROUTINE TMP(ND) 

IMPLICIT REAL*:=8( A-H,Q-Z) 

C0MMCN/8L1/P< 51,51) ,G(51,51),SI(51,51) 

CCMHON/BL2/BFT PHI , B Z , NCOL S , NROWS , NCE N , MA X 
C0MMCN/BL6/TH(51 ,51 ) ,DTX(51 ) 

NRCM=NROWS-l 
NCCM=NCCLS-1 
G1 =GAM>:<.25D 00/(H-H) 

G2=GAM/ (H*H ) 

DO 2 I=2,NRCM 
DO 2 J=2,NCCM 
EL = G( I , J) 

DPX=P { I + l , J )-P ( I-l , J ) 

DPY=P(I , JM )-P( I , J-1 ) 

DNX=G( I+l , J )-G( I-l , J I 
DNY = G ( I , J + 1 )-G( I , J-1 ) 

ARG = G1*( ONX^‘DPX+ONY-DPY)/EL-l.D 00 
ARG2 = G2-( DPX^-DPXfOPY=:=DPY) 

TH(I , J}=< -ARG + DSQRT( ARG=^ARG + ARG2) )/2 .D 00 

2 CONTINUE 
IF(NC.EQ.l) GO TO 5 

TF = (1.D 00 + DSQRT(l.D 00+4«D 00 «GAM# P ( NROM, NC EN ) -*2 ) ) / 2 
.D 00 

DO 3 J=l, NCOLS 
THINROWS, J ) =TF 

3 Th(l,J) = UD 00 
DO 4 I=2,NR0M 

Trid ,NCOLS)=TH(I ,1) 

4 TH(I,1) = TH( I,2)-THU,NC0LS)+TH(ItNC0M) 

RETURN 

5 TF=(1.D C0+DSQRT(1.D 00+4, D 00«G AM* P ( NCE N , N COM ) *«2 ) ) / 2 

,D 00 

DO 6 I=1,NRCWS 
TH( I , 1 )=1 ,D 00 

6 TH(I , NCOLS! =TF 
DO 7 J=2,NCCM 

TH (NRCN'S , J) =TH( 1 , J) 

7 TH(1 , J)=TH( 2, J)-TH(NROKS, J) +TH(NROM, J) 

RETURN 

END 
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CURRENT PROGRAM 






c 



c 



c 



THIS PPCGPAM IS DESIGNED TO CALCULATE THE CURRENT 
DISTRIBUTION GIVEN THE POTENTIAL AND CHARGE DENSITY 
FIELDS. IT OUTPUTS USING SUBROUTINE WRITER FROM TFE 
SHEATH PROGRAM. AS WITH THE SHEATH PROGRAM IT USES A 
DATA CELL FOR VARIABLE STORAGE. IT REQUIRES THE 
FCLLCKING PREDETERMINED PARAMETERS: 

P POTENTIAL DISTRIBUTION 

G ELECTRON DENSITY DISTRIBUTION 

SI ION DENSITY DISTRIBUTION 

GAM JOULE HEATING PARAMETER 

BET HALL PARAMETER 

ZL STEP SIZE MULTIPLE 

NROWS ARRAY WIDTH 

NCOLS ARRAY LENGTH 

IMPLICIT REAL=.‘'3( A-H, 0-Z ) 

DIMENSION P(5l,5l),G(51,5l),SI(51,51),XJ<51,51),YJ(5l, 
51), SX(51,51»,ZI(51),TH (51,511 
READ(l) P 
READ(l) G 
READ(l) SI 
GAM=0.0 
ZL=5.D 02 
BET=1.D 00 
NRGWS=51 
NCCLS=51 
NRCM = iNRCWS-I 
NCCM=i\CCLS- I 
Nd=NCCLS-l 



NC=NB-1 

H=ZL/DFLCAT ( NROWS- 1 ) 

HR=5.0C“01/H 
G1=GA(->^.?5D 00/(Fp:^H) 

G2=GAf-'/(H’:'H ) 

DO 10 1=1, NROWS 
DC 10 J=l, NCOLS 
10 Th( I , J ) =1 .0 00 

IFCGAM.LE.l .D-04) GO TO 20 
CALCULATE DERIVATIVES AND ELECTRON TEMPERATURE 
DC 2 I=2,NR0M 
DC 2 J=2,NCCK 
EL=G( I , J) 

DPX= F( H-1 , J )-P( I-l , J ) 

DPY=P ( I , J^l )-P( I , J-1 1 
DNX=G(I+1 ,J)-G( I-l, J ) 

DNY=G( I , J+L )-G( I , J-l ) 

ARG=Gl-( DNX'^DPX+-DNY=i^DPY )/ EL-1 . 0 00 
ARG2 = G2^'( 0FX=i=0PX4-DPY-DPY) 

TH(I , J) = (-ARGi-DSQRT( ARG*ARG+ARG2 ) )/2.D 00 

CALCULATE^CROSS-CHANNEL CURRENT DENSITY 

SI ( I , 1 )=-3!^*Gf i,l)+4.*G(I,2)-G(I,3)-G(I,U*(- 
3.>?P( I , 1 )+4.-P( I , 2)-P( 1,3) ) 

SI (I , NCOLS )=G( !,NC)-4.-G( I,NBM-3.*G( I, NCOLS )- 
G(I ,NCOLS)*( P( I,NC)-4.«P( I ,NB) + 

3.«P( I , NCOLS) ) 

3 si(i, J)=Tm I, j)*(G(i , j+i)-G(i, j-1) )-G(i, J)-{ p(i ,J<-n- 

CALCULATE^AXIAL^CURRENT DENSITY 

SX(l,J)=-3^^Gfl,jM-4.*G(2,J)-G(3,J)-G(l,J)*{- 

3.’^P( l,J)■^4.*P(2,J)-P(3,Ji ) 

SX(NRCWS, J)=G(NC, J )-4.«G(NB, J) <-3.=^G(NR0WS, J )- 
G (NROWS , J )*( P(NC , J )-4 .*P( NB, J ) <■ 

3.*P( NROWS, J) ) 

4 SX(t,J)=TH( I,J)=!'(G( I+1,J )-G(I-l,J) )-G( I,J)*( P( I + 1,J)“ 

P( I-l , J) ) 
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C APPLY HALL PARAMETER 
DC 5 I=1,NR0WS 
DO 5 J=1,NCCLS 

YJ(ItJ) = (SI(I,J)-BET«SX(I,J) )’^HR 
XJ(I , J)=( SX( I, J)+BET*SI( I , J) )*HR 
5 CONTINUE 

C INTEGRATE FOR NET CURRENT 

FAC=1.0 00/(3. D 00=^'DFL0AT( N6) ) 

DO 7 I=1,NR0WS 
SUPl =C .0 
SUM2=C.O 

SUM=XJ(I, 1 ) +XJ(I ,NCOLS> 

DO 8 J=2,NB,2 

8 SUM1 = SUM1<-X J( I , J) 

DO 9 J=3,NC,2 

9 SUM2=SUM2+X J( I , J) 

ZI (I ) = ( SUM+4.D 00>'-SUMl + 2. 0 00=(< SUM2 ) 4 F AC 
7 CONTINUE 
C OUTPUT 

WRITE(6, 101 ) 

CALL WRITER(X J , NRCK'S ,NCOLS ) 

WRITE(6,i02) 

CALL WRITER (YJ,NROWS,NCOLS) 

WRITE (6, 104 I (1,21(1), 1 = 1 ,NROWS ) 

STCP 

101 FORMAT (//,3X, 'X-CURR5NT* ,/) 

102 FORMAT (//, 3 X, ' Y- CURRENT* ,/ ) 

104 F0RMAT(//,51( 3X, 14, 1P012.5,/) ) 

END 
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c** 



CNE-DIMENSIONAL POTENTIAL AND CHARGE DISPLAY 






C 



C 



c 



THIS PROGRAM DTSPLAYS THE POTENTIAL AND CHARGE DENSITY 
DISTRIBUTIONS ALONG A CUT PERPENDICULAR TO THE WALL 
AND THROUGH THE ELECTRODE. IT DISPLAYS THE RESULT AS 
GRAPH ON THE CALCOMP PLOTTER USING A NAVAL POSTGRAOUAT 
SCHOOL SLBROUTINE ORAWP. IT REQUIRES THE FOLLOWING 
PREDETERMINED PARAMETERS: 

P POTENTIAL DISTRIBUTION 

G ELECTRON DENSITY DISTRIBUTION 

SI ION DENSITY DISTRIBUTION 

RAT CHARGE DENSITY BOUNDARY VALUE 

NM ARRAY LENGTH 

NC ELECTRODE LOCATION 

ZL STEP SIZE MULTIPLE 

TITLE DESIRED GRAPH LEGEND 

OTHER PARAMETERS ARE EXPLAINED IN SUBROUTINE DRAWP. 
IMPLICIT REALMS (A-H,0-Z) 

INTEGER^^'A I T 3 U 2 ) / 1 2 ^J=0/ 

P,EAL*4 RTB ( 28 I /Za^^O. 0/ 

REAL*A LABL(3}/«P0T ‘,‘POS ',’NEG »/ 

DIMENSION X(5U ,Y(51 ) ,Z(5n,VAl(5U 
DIMENSION A(51t51) 

EQUIVALENCE ( T I T LE , RT B ( 5) ) 

REAL*8 TITLE! 12) 

READ(5,101) TITLE 
SET CCNSTANTS 
RAT=1 .0-03 



NM = 51 
NC=26 
ZL=5.D C2 
H=ZL/CFL0AT(NM-1 ) 

ITB(4)=10 
RTB{ 2)=. 125 
I TB( 5)=l 
IT6( 6 1=1 
ITB( 8)=2 
DO 12 1=1, NM 
12 X( I)=H=:'DFL0ATU-1) 

INPUT DATA 
READ(l) A 
AA=A(1 ,26 ) 

DC 1 1=1, NM 

1 Y ( n=A( I ,NC )/AA 
READ(l) A 

DC 2 1=1, NM 

2 Z( I)=A(I,NC) 

REAOm A 

DO 3 1=1, NM 

3 W(I)=A(I,NC) 

NFLAG=1 

10 GO TO (4,5,6,7) ,NFLAG 

5 DO 8 1=1, NM 

8 Y(I ) =W( I > /RAT 
GO TO 4 

6 DC 9 1=1, NM 

9 Y( I)=Z(I)/RAT 

4 RTB( 3) =LABL{NFLAG) 

ITB( 1 )=MFLAG 

OUTPUT 

CALL CRAWP ( NM, X, Y ,I TB ,RTB) 
NFLAG = NFLAG*-1 
GO TO 10 

7 STCP 

101 FCPMAT(6A8) 

END 
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ELECTRON TEMPERATURE PROFILE 



* 4 : 



C 



C 



C 



THIS PROGRAM CALCULATES THE ELECTRON TEMPERATURE RATIO 
ALONG A CUT PERPENDICULAR TO THE WALL AND THROUGH THE 
ELECTRODE. IT DISPLAYS THE RESULT AS A GRAPH ON THE 
CALCGMP PLOTTER USING THE NAVAL POSTGRADUATE SCHOOL 
SUBRCUTINE DRAWP. IT REQUIRES THE FOLLCWING 
P R E - D E T E R M I N E D PARAMETERS; 

P POTENTIAL DISTRIBUTION 

G ELECTRON DEI'ISITY DISTRIBUTION 

SI ION density DISTRIBUTION 
GAM JOULE HEATING PARAMETER 
Z STEP SIZE MULTIPLE 

NH ARRAY LENGTH 
NC ELECTRODE LOCATION 
TITLE DESIRED LEGEND GN GRAPH 
OTHER PARAMETERS ARE EXPLAINED IN SUBROUTINE ORAWP 
IMPLICIT REALMS {A-H,0-Z) 

INTEGEP- A I TB( 12J/12=i=0/ 

REAL-"4 PTB( 28i/23=<=0.0/ 

PI PENSION P( SI ,51 ) ,G( 51,51 ) ,SI ( 51, 51 ) ,TH(51 ) ,XI 51 J 
EQUIVALENCE ( T I T Li: , R T 8 ( 5 ) ) 

REAL^-8 TITLE! 12) 

READ(5,101) TITLE 
SET CONSTANTS 
GAM=1.0D 03 
Z=5.CC 02 
NM = 51 
NRCM=50 

H=Z/DFLQAT(NROM) 

MC = 26 
NCM=NC-1 



K'CF=NC + 1 

G1=GAM=<--.25D 00/(H=:=H) 
G2=GAM/(H«H) 



IT B ( 4 ) = 1 0 
ITB( 5 ) = i 
I TB( 6)=1 
ITB( 8 ) = 2 
RFAD(l) P 
READ(l) G 
PEAOm SI 
DO 12 1=1, NM 

12 X( I ) =H^^DFL0AT( I- 1) 

TH( 1 ) =1 .D GO 
TH(NM)=1.D 00 

CALCULATE ELECTRON TEMPERATURE RATIO 
DO 1 I=2,NR0M 
EL=G( I ,NC) 

DPX=P ( I +1 ,NC) -P( I-l , NC) 

DPY=P (I,NCP )-P(I ,NCM) 

DNX=G( 1+1 ,NC) -G< I-l ,MC) 

DNY=G(I ,NCP)-G( I ,NCM) 
ARG=G1=<^{DNX':=DPX+DNY"0PY)/EL-1.D 00 
ARG2 = G2*( DPX=’?DPX + CPY*DPY) 

THI I ) = (-ARG+OSCRT( ARG-ARG+ARG2) )/2.D 



00 



1 CONTINUE 
OUTPUT RESULTS 

WKITE(6,iC0) (THU) , 1 = 1 ,NM) 
CALL CRAWP(NM,X, TH, ITBtRTB) 
STOP 

100 F0PMAT(51I3X,F12.3,/ ) ) 

101 FCRMAT(6A8I 
END 
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CONTOUR PROGRAM 






THIS PROGRAM PLOTS THE CONTOURS OF THE POTENTIAL AND 
CHARGE DISTRIBUTIONS USING THE CALCGMP PLOTTER. IT USES 
THE SAME DATA CELL AS THE SHEATH PROGRAM. THE 
FOLLOWING PARAMETERS MUST BE PREDETERMINED: 

AD POTENTIAL OR CHARGE DENSITY ARRAY 

NRCWS ARRAY LENGTH 
N'CCLS ARRAY WIDTH 

NL NUMBER OF POTENTIAL CONTOURS TO BE PLOTTED < ZI 

MLl NUMBER OF CHARGE DENSITY CONTOURS TO BE PLOTTED 
IW CESIRED OUTPUT WIDTH IN INCHES < 10 
IH DESIRED OUTPUT HEIGHT IN INCHES < 16 
TABL 6 FORMATTED VARIABLES 

(1) JOULE HEATING PARAMETER 

(2) TEMPERATURE 

(3) ELECTRODE POTENTIAL 

(4) CHARGE DENSITY BOUNDARY CONDITION 

(5) HALL PARAMETER 

(6) CHARACTERISTIC LENGTH 

LTG 3 GRID SELECTIONS TRUE/FALSE 

(1) TRUE. ..ALL EXTERIOR CONTOUR SEGMENTS AND 

INTERNAL MAX I MUMS WILL BE LABELED. 

(2) TRUE ... REQUEST TIC MARKS AND CORRESPONDING 

SCALE VALUES ONE INCH APART ON EXTERIOR 
FRAME OF THE CONTOUR GRAPH 
FALSE.. OMIT TIC OPTION 

(3) TRUE. . .REQUEST A ONE INCH BY ONE INCH STRAIGHT 

LINE GRID SUPERIMPOSED ON THE CONTOUR 
GRAPH 

FALSE.. CM IT GRID 

THE FROGFAM HAS THE FOLLOWING SUBROUTINES: 

MAIN. . .GENERAL COORDINATION AND PREPARATION OF 
OUTPUT 

FLOP. ..ORIENTS ARRAYS TO CORRESPOND WITH OUTPUT 
LABELLING CONVENTION 

CLVL. . .DETERMINES LINES OF CONSTANT VALUE TO BE 



CCNTUR, 



PLOTTED 

..LABELLING AND PLOTTING COORDINATION, 
ALLOWS SUPERPOSITION OF 3 CONTOURS ON ONE 
PLOT. IT UTILIZES ADDITIONAL SUBROUTINES 
SUCH AS RESTOF, SCAM, TRACE, CALC AND PLOTT 
AND IS A MODIFICATION OF AN NPS PROGRAM. 
CGMMCN/BL 6 /NRO WS , NOG L S , I W , I H , N 1 
COMMCN/BLl 1 /XLIN( 100) , YLI N( 103) 

REAL=i^8 TITLE(30)/' • , * ' » 

I • 

•DENSITY • 

I I 

( I 

•LENGTH • 

• ON • , 

•BOX 1219^/ 

• CHARGE 



1 • 

2^ 

5 • 

4^ VOLTAGE 

5 • 

6 ' 

?• AL MAP 
REALMS 



TEMP 



R. DDLS 



, ■ GAM 

,• 

,• 

, •BET 

,• 

• POTENTI 



TITLFA(4) /• ELECTRON^ , 
REAL'^8 TABLE(6) 

REAL^E AD(51,51) 

LOGICAL*! LTG(3) 

DIMENSION A(81,61) 

DIMENSION CL(30) 

NAMELI ST/TA6L/TABLE , LTG 
NAMEL IST/OOLL/NROWS,NCOLS 
READ ( 5, TABL ) 

READ( 5,D0LL ) 

XVAL=(IW*1. )/(NCOLS*l.-l. ) 
YVAL=( IH*1 . )/(NROWS*l .“1. ) 



• , • 



ION CHAR' , • GE 



NL,NL1, IW, IH 



N1=0 

DO 10 1=1,1 
N1 = NH-1 

XLIN( M)=XVAL*(20) 
10 YLIN(M)=0.0 
DO I K=l,3 
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!’ 

i 






r 






C INPUT ARR/iY 

READI2) AO 
DO 19 I=l,NROWS 
DO 19 J=l,NCOLS 
19 AU, J)=AO( ItJ) 

CALL FLCPIA ,NROWStNCOLS) 

C FIRST FIGURE CONTOUR LEVELS TO BE PLOTTED FOR 

11 CALL CLVL( A,CL,NL) 

WRITE(6,210) 

DO 3 I=1,NL 

3 WRITE(6,215) I ,CL( I ) 

12 TI 7LE(8)=TABLEm 
TITLEUl > =TABLE(4) 

TITLE! 14)=TABLE( 2) 

TITLE!17)=TABLE(3) 

TITLE(20)=TA3LE( 5) 

TITLE! 23) =TABLE! 6) 

CALL CCNTUR ! A jNROWS,NCOLS,81 ,CL,NL ,TITLE» LTG ) 
C NOW FLCT CHARGE FIELD CONTOURS 

IF!K.EQ.3) STOP 
KL = 2=!^K-1 
KF=KL+1 

TITLE(27) =T ITLEA!KL) 

TITLE(28> =TITLEA!KP) 

1 CONTINUE 

210 FGPMAT!' «,T26,'THE LEVELS PLOTTED*) 

215 FORMAT!' • , T33 , * CL ! ' , I 2 , ‘ ) =* , F8 .4 ) 

STOP 

ENC 



SUBROUTINE FLOP! Z , NR CWS , NCOLS ) 
D I PENS ION Z(8l»6l) 

I I VRT=NRuWS/2 
DO 3 1=1,1 IVRT 
M=NR0NS-!I-1) 

DO 3 J=1,NCCLS 
SAVE=Z! I , J) 

Z ! I, J) = Z! M, J ) 

3 Z!P,J)=SAVE 
RETURN 
END 



SUBROUTINE CLVL ( CM, C LM , NU ML ) 
COPMCN/BL6/NROV,S , NCOLS , IW , IH,N1 
DIMENSICN CM! 81 ,61 ) ,CLM( NUML) 

CMIN=CM! 1,1) 

CMAX=CP IN 
DO 5 J=l, NCOLS 
DO 5 I=1 ,NkCWS 

I F !CMU , J ) . LT oCM IN ) CM IN=CK! I , J ) 

IF!CM! I , j; .GT.CMAX) CHAX=CM! I , J ) 

5 CONTINUE 

C NON DETERMINE CONTOUR LEVELS TO BE PLOTTED. 

J=NUML 
I=NUML-1 
ANL= 1=^-1 . 

PLIMT=!CMAX-CMIN)/ANL 
C NOW FILL THE CONTOUR LEVEL VECTOR 

DO 6 1=1, J 

6 CLP! I )=CHIN«-! I-l ) -'PL INT 
RETURN 

END 



SUBROUTINE CONTUR! AM ,M,N, MX ,CL, NL, T ITLE,LTG) 



COMMON /INTFAC/ X,Y 
COMMON /DAYFOF/ MT,NT,NI 



IX, lY , IDX, IDY, ISS, IT, 



MATRIX A 



IV,NP, 
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1 

60 

71 

64 

2 

40 

61 

3 

4 

5 



C0yMGN/TA3L/ T ABC ( 20 , 6 ) , JC 

COMMCN/OITS/XMIN,YiMIN, SLOPEX, SL0PEY,DITS0X, DITSDY, 
C0MMCN/BL6/NR0vvS ,NCOLS , IW, IH,N1 
C0MY!CN/BL11/XLIN(100 1 ,YLIN( 100) 

CGMMC.\/8L13/\'SCN 
REAL=’^8 TITLE(l) 

WIDTH/' WIDTH* /, HE IGHT/ 'HE IGHT' /, WHICH 
DIMENSION AMI MX, 1) »CL ( 1) 

DIMENSION RECI900), XI1800), YI1800) 

DIMENSION IPT(3,3) , INXI8) , INYI8) 

DIMENSICN OITSXI 5 ) , 01 TSYI 5 ) 

LCGICAL^n LTGI 1 ) ,MINUS,LABL 
JC = 0 

LABL=LTG(1 ) 

CHECK IW PARAMETER 
WHICH=YIDTH 
IFIIW) 1,1,2 
WRITE(6,50) WHICH 

FORMAT! 'O' ,T7,A8, 'OF CONTOUR GRAPH ILLEGAL.') 

WRITE (6,54) 

FORMAT! 'O’ , T7, 'NO GRAPH WILL BE PRODUCED.' I 
RETURN 

CHECK IF IW IS TOO WIDE 
IF(IW-9) 3,3,40 
WRITE! 6, 61) 

FORMAT! 'O' , T7, ' I W PARAMETER GREATER THAN 9. CONTUR 
I W = 9 

NOW CHECK IH PARAMETER 
I F! IH) 4,4,5 
KHICH=HEIGHT 
GO TO 1 

DITSDX=!N-1 ,0)/IW 
0ITSDY = !-1 .0^-M)/IH 



XMIN=1 .0 
YMIN=-M 

SLCPEX=1.0/DITSDX 
SLGPEY=1.0/CI TSDY 
DI TSX! 1 ) = 1.0 
0ITSX!4)=1.0 
DITSX! 5)=1 .C 
DITSXI2 )=N 
DITSX!3)=N 
OITSY! 1 )=-l .0 
DITSY!2)=-1 .0 
DITSY!5)=-1.0 
0ITSY!3)=-M 
0I7SY!4)=-M 
DC 201 1 1=1,5 

DITSX! I)=SL0PEX*!0ITSX! I )-XMIN) 

DI TSY! I ) = SL0PEY* IDITSY! I )-YMIN) 

STARTP=(9.0-IW)/2.0 

CALL RLOTS 

CALL PL0T!STARTP,0.0,-3) 

CALL LINE!DITSX,DrrSY,5,l,l) 
DITSX! 1 )=DITSX( 1 )-.5 
DITSX!5)=DITSX(1 ) 
niTSX!4)=DITSX(4)-.5 
DITSX! 2)=DITSX!2)+.5 
DI TSX(3)=DITSX!3 )+-.5 
DITSY! I)=OITSY! 1 )+.5 
DITSY!5)=DITSY! 1) 

DITSY!2)=DITSY!2 )+.5 

DITSY!3)=DITSY(3)-.5 

0ITSY!4)=DITSY(‘t)-.5 

CALL LINE!0ITSX,DITSY,5,1, 1) 

SLCPEX=1 .0/DI TSDX 

SLCPEY=1.0/CITSDY 

I ENDX=SL0PEX-N*1 

I ENDY=SLCPE Y«M•^1 

IF! .N0T.LTGI2 ) ) GO TO 34 

DRAW TIC MARKS ON OUTtR FRAME 

START CN LFFT EDGE GOING DOWNWARD 
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IFLAG=0 
ZINGX=-. 1 
Z IKGV=C.O 
ZX=0.0 
Z Y ~ ^ 

cx=nffsx( 1 ) 

CY = 0 ITSY( L )-.5 
I £ND= 1 1 

2222 IFLAG= IFLAG + 1 

DC 2022 1=1 ,I END 
CALL PLOT (CX,CY.3) 

COCRDX=CX+Z INGX 
CCCkDY = CY4-Z INGY 
CALL PL0T(CG0RDX,C00RDY,2) 

CX=CXf ZX 
2022 CY=CY+ZY 

GO TO <21»22,23,24) , IFLAG 
C NOW DO THE RIGHT EDGE GOING DOWNWARD 

21 ZINGX=.l 

CX = DIT SX(2 I 
CY = OnSY( 2)-.5 
GO TO 2222 

C NCW DC TOP EDGE 

22 ZINGX=C.O 
ZINGY=.l 
ZX=. 9 
ZY=O.C 

CX=D ITSX( 1) f .5 
CY=D ITSYI 1 ) 

IEND=11 
GO TC 2222 

C NOW DC THE BOTTOM EDGE 

23 ZINGY=-.l 

CX = OITSX( 4) +.5 
CY=0ITSY(4 ) 

ZINGY=-. 1 

GC TC 2222 

NOW LABEL TIC MARKS 

DO X-CiRECTION FIRST, TOP EDGE 

POSITION PEN 

24 DELTAX=.l 
IFLAG=C 
ZX = . 9 
ZY=O.C 

CX = DHSX( 1 )+.35 
CY = DITSY( U + .12 
3033 IFLAG=IFLAG+1 
XZERC=0.0 
DO 3333 1=1,11 

CALL NUMBER (CX,CY, . 1 4 , XZE RO ,0 . 0 , 2 ) 

CX=CX+ZX 

CY=CY+ZY 

3333 XZERO=XZERO+DcLTAX 

GO TO (31 ,32,33,34) , IFLAG 
C LABEL BOTTOM EDGE TIC MARKS 

31 CX=DITSX(4H-.35 
CY = 0 ITSY(4 )-.19 
GC TO 3033 

C LABEL LEFT EDGE OF TIC MARKS 

32 CX=D ITSX(4)-.4 
CY=D ITSY( 41+.46 
DELTAX=. 1 

I END= I ENDY 
ZX=0.0 
ZY = . 9 

GO TO 3033 

C NOW LABEL RIGHT EDGE TIC MARKS 

33 CX=D IT5X( 3 ) + .12 
CY=DITSY(3)+.46 
GC TO 3033 

C CHECK IF GRID DESIRED 

34 CALL RESTOF(LTG, lENDX, IENDY,NL, AM,M,N,MX,CL, STARTPtTIT 
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c 

c 



4044 



4 444 



42 



35 

20 



777 

778 



779 



LEiOITSX, DITSY) 

RETURN 

END 

SUBROUTINE RES TOP ( LTG , I ENOX , I ENOY t NL , AM t M, N , MX » CL , ST AR 
TPtTITLE, DITSX, DITSY) 
realms TITLEd) 

DIMENSION AM(MX,1) ,CL( 1) 

DIMENSION QI TSX( 5) ,01 TSY( 5) 

LOGICAL-l LTG( 1) , MINUS, LABL 
C0MMCN/TABL/TABC(20,6) , JC 

COVMON/DI TS/XMIN,YMI N, SLOFEX, SLOPEY, DITSDX, DITSOY, 
C0MM0N/BL6/NR0WS,NCQLS ,I W, IH,Nl 
COMMCN/BLll/XLIN(iOO) , YLIN( 100) 

C0MMCN/BL13/NSCN 

IF( .NCT.LTGO) ) GO TO 35 

DRAW INCH BY INCH GRID 

I ENO= IENDX-2 

POSITION PEN 

IFLAG=C 

CX=D ITSX( 1 ) + .5 
CY = DITSYa I-.5 
COCROX^O.O 
CCGRDY=-IH 
DX = 1 .0 



DY=0.0 

DO 4444 1 = 1 , I END 
CX=CX+DX 
C Y=f Y + r Y 

CALL PLOT (CX, CY,3) 

ZX=CXtCOCRDX 

ZY=CY+CGOROY 

CALL PLCTIZX, ZY,2J 

IFUFLAG) 35,42,35 

IFLAG=1 

IEND= IENDY-2 

CY=DITSY(4H-.5 

CX=DITSX(4)+.5 

CCCRDX=1W 

CCCRDY=0.0 



DX=0.0 
DY = 1 .0 
GO TO 4044 
CONTINUE 

CALL L I NE ( XL I N , YL I N , N 1 , 1 , -6 ) 

DO 20 1=1, NL 

CALL SCAN(AM,M,N,MX,CL( I) ) 

NSCN=1 

IF(. NOT. LABL) GO TO 778 
IF(JC.EQ.O) GO TO 778 
DO 777 1=1, JC 
CGCRDX=TA3C( 1,4) 

COCRDY=TABC{ 1,5) 

CLEV=TA8C( I ,6 ) 

CALL NUMBER (COORDX,COQRDY , . 07 , C LE V , 0 . 0 , 3 ) 

CONTINUE 

CALL SYMBOL { -ST ART P, IH^-1.0,.21,TITLE(25),0.0,46 ) 

CALL SYMBOL ( -5TARTP , IH + 1. 5, .21, TITLE! 19), 0.0, 48) 

CALL SYMBOL (-STARTP, IH+-2.0 , . 2 1 , T I TL E ( 13 ) ,0 . 0 , 48 ) 

CALL SYMBOL (- ST ARTP ,IH+2.5,.21,TITLE(7).0.0,48) 

CALL SYMBOL! -STARTP, IH+-3.0,. 21, TITLEd), 0.0, 48) 

CALL FLCT(-STARTP,IH+6.5,-3) 

CALL PLOTE 
RETURN 

CALL SYMBOL ( -ST ARTP, IH+1.0,.21,T[TLE(25),0.0,48) 

CALL PLOT(-STARTP, IH+4.5,-3) 

CALL FLCTE 

RETURN 

END 

SUBROUTINE SCAN! AM, M , N ,MX , CL ) 

DIMENSION AM(MX,1 ),REC(900), X!1300), Y(1800) 

DIMENSION IPT!3,3) ,INX(8) , INY!8) 

COMMON /OAYHOF/ MT , NT , N i , I X , I Y , I DX , I DY , I SS , I T , I V , NP , N G 
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PY, REC,CV, 



I FT, INX, 



, J T , 

INY,DL,RA,THE 
COMMON /I NT FAC/ X,Y 
LOGICAL -n LA8L, MINUS 
COMMON/OITS/XMIN,YMIN,5LOPEX, SLOPEY , D ITSDX , D IT S CY 
COMMON/6L 13/NSCN 
D = 0. 

R = l. 

TH = 1.57C796 

NP=0 

DL=D 

RA=R 

TFE=TH 

MT = N 



NT=M 

CV=CL 

I FINSCN.EQ. 1 ) IZW=0 
I F( I ZW'120631 ) 1,3,1 

1 IPT(1,1)=8 
IPT( 1 , 2)=1 
IPT( 1 ,3)=2 
I PT(2,1)=7 
IPT( 2,5)=3 
I PTI 3, 1 )=6 
I PT( 3, 2)^5 
IPT(3,3)=4 
INX( 1 )=-l 
INXI 2) =-l 
INX( 3)=0 
INX(4)=1 
INXI 5 )=1 
INX(6)=1 
INXI 7)=0 
INXI 8 )=-l 

I NY I 1 )=0 
INYI 2 )=1 
INYI 3»=fl 
I NYI 4 )=+l 
INYI 5 )=0 
INYI 6)=-l 
INYI 7 )=“1 
INYI 8)=-l 
I ZU=120631 
3 XT=MT 

DO 58 J=l,900 
58 RECI J )=0 
I 5S = 0 

2 MT1=MT-1 
IDIR=1 

DC 110 1=1, MTl 
IFIAMI 1,! I-CV) 55,110,110 
55 IF I AM I 1, I +-1 )-CV) 110,57,57 
57 IX = H-1 



IY=1 

IDX=-1 

I0Y=0 

CALL TRACE I AM, MX) 

110 CONTINUE 
NT1=NT-1 
I DIR=2 

DC 20 1=1, NTl 
I FIAMI I,MT)-CV) 15,20,20 
15 IFIAMI I+1,MT)-CV) 20,17,17 
17 1X=MT 
IY=I +1 
IDX=0 



IDY=-1 

CALL TRACE I AM, MX) 
20 CONTINUE 
IDIR=3 

22 DO 30 1=1, MTl 
MT2 = MT^-1-I 



, 
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IF( AM(M,HT2)-CV) 25,30,30 
25 IF(A/^<NT,MT2-1)-CV) j0,27,27 
27 IX=MT2-1 
I Y=NT 
IDX=1 
IDV=0 

CALL TRACE (AM, MX) 

30 CGNTIKLE 
IDIR=A 

DO 40 1=1, NTl 
M2 = M + 1 -I 

IF( AM (NT2, n-CV) 35,40,40 
35 IF(AM(NT2-1 , 1 )-CV) 40,37,37 
37 IX=1 

IY=NT2-1 

I0X=0 

IDY=1 

CALL TRACE (AM, MX) 

40 CONTINUE 
IDIR=5 
ISS=1 
N'T 1 = NT-1 
MT1=MT-1 
00 10 J=2,NT1 
DO 10 1=1, MTi 
IF( AM{ J, I )-CV) 5, 10, 10 
5 IF (AM( J, n-1 )-CV) 10,7,7 
7 C0M=1C0=<^( I + l ) i-J 
IF (NPJ 12,11,12 
12 00 9 ID= 1, NP 

IF (RECUD)-COM) 9,10,9 
9 CONTINUE 
11 IX= I+l 
IY = J 
IDX=-1 
IDY=0 

CALL TRACE (AM, MX) 

10 CONTINUE 
RETURN 
END 

SUBROUTINE TRACE (AM, MY) 

DIMENSION AM(MY,1 ), RFC (900), X(1800), Y(1800) 

D I ME NS I CN IPT(3,3), I N X ( S ) , I NY ( 8 ) 

COMMCA' /DAYHOF/ MT , NT , N I , I X , I Y , I OX , I DY , I SS , I T , I V , NP , 

COMMON /INTFAC/ X,Y 

PY=0.0 

RC= CCS (THE)-RA 
RS= SIN (THE)=^RA 
501 JT=0 
N = 0 

IXC=IX 

IYO= lY 

ISX= I OX+2 

ISY=I0Y+2 

IS=IPT( ISX, ISY) 

JTB = 0 
ISO=IS 

IF(ISC-8) 18,18,17 

17 ISC=ISC-8 

18 IT=0 

5 CONTINUE 

CALL CALC (AM, MY) 

NZ = N 
N = N2 

IF ( IT + JT-1 ) 49,49,47 
47 XS=X(N-1) 

YS=Y(N-1 ) 

X( N-1 )=X(N!) 

Y(N'1)=Y(N) 

X( N) =XS 
Y(N)=YS 
49 IS=IS+1 
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I 



9 

7 

8 

308 

103 

51 

20 

21 

22 

23 

10 

13 

19 

11 

12 

206 

213 

217 

214 

215 

6 

306 

16 

50 

307 

73 

74 

1 

2 

41 



JT=1T 

IF ( IS-9) 8,7,7 
I S=I S-8 
IDX=IKX( IS) 

IDY= INY( IS ) 

1X2= I X+IDX 
IY2= I Y+IDY 
JT6= JTB+1 

IF {JTB-1799) 51,51,308 

PRINT 103,CV,X( N) ,Y(N) 

FGPMAT(1H0,23HA CONTOUR LINE AT lEVEL,E12.5, 

RETURN 

CONTINUE 

IF ( I SS) 10,10,20 
IFMX-IXO) 12,21,12 
IF(IY-IYO) 12,22,12 
IF(IS-ISO) 12,23,12 
CONTINUE 

CALL CALC (AM, MY) 

GG TO 73 

IF(IX2) 13,50,13 

IF (IX2-MT) 19,19,50 

IF ( IY2) 11,50,11 

IF (IY2-NT) 12,12,50 

I F(CV-AV|( IY2, 1X2) ) 206,206,5 

IF ( I CX=:=*2+ IDY^‘==;^2- 1 ) 213,6,213 

DCF= ( AM( lY, IXH-AM( lY , 1X2 ) + AM( IY2, IX ) fAM( IY2 , 1X2 ) ) /4.0 

IF (DCF-CV) 5,217,217 

IF (INX(IS-ll) 214,215,214 

I X=I X+IDX 

IDX = - IDX 

PY=2 . 0 

CALL CALC (AM,HY) 

IX=IX+I0X 
GG TO 6 
IY=I Y+IDY 
IDY=“ IDY 
PY=2.C 

CALL CALC (AM, MY) 

IY=1 Y+IDY 

I F( AMI lY, IX-1 )-CV) 306,15,16 
NP-NF+ I 

REC( NP) = 100=!=IX + IY 

I S=I S+5 

1X=I X2 

IY=IY2 

GO TO 9 

XT=MT 

I F(AM( lY, IX-1 )-CV) 307,73,73 



NP=NP+ 1 

R.EC( NP) = 1C0'-:^IX + 1 Y 

DO 74 1=1 ,N 

X( I ) = X( I ) + RC*Y(I ) 

Y( I )=RS’;^Y( I ) 

CALL PLCTT(N,CV) 

RETURN 

END 

SUBROUTINE CALC(AM,MY) 

DIMENSION AM(MY,1 ),REC(900), XdSOO), Y(1800) 
DIMENSION IPT(3,3) ,INX(8) ,INY(8) 

COMMON /OAYHOF/ MT , N T , N I , IX , I Y , I DX , I OY , I SS , I T , I V , NR , N 
COMMON /INTFAC/ X,Y 
I T=0 



N = N+1 

IF ( IDX*«2 + IDY=!=«2 -1) 20,1,20 
IF (IDX) 10,2,10 
X(N) =IX 
Z= lY 

IY2= lY+IDY 
DY=IDY 

Y(N) = ((AM(IY,IX)-CV)/(AM.( I Y , I X ) - AM ( I Y 2 , I X ) ) )«DY+Z 



RETURN 



, 
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c 



c 



c 



c 



10 



44 

20 



24 

21 

23 

27 



25 

33 

28 



100 



1 

2 

5 

3 



4 



6 



200 



Y(N)=IY 
W= IX 
DX=IDX 
1X2= IX+IDX 

X (K) = ( ( AM( lY, IX)-CV)/{ AM( lY, IX) -AMI I Y, 1X2) ) ) *DXtW 

RETURN 

1X2= IX+IDX 

IY2= I V+IDY 

W= IX 

Z=IY 

DX=IDX 

DY=IDY 

DCP=( AMI lY, IX) 4-AMI IY t 1X2) +AMI IY2, IX ) + AMI IY2 f 1X2 ) ) /4.0 

IF IPY-2.0) 24,21,24 

IF IDCP-CV) 21,21,25 

AL=AMI lY, IX)-DCP 

V= .5«I AL^-DCP -CV)/AL 

XIN) =V>^-DX<-W 

Y I N)=V*0Y+Z 

PY=0.0 

RETURN 

IT = 1 

AL=AMI IY2, I X2 ) -OOP 
V=,5=i=IAL + DCP-CV) /AL 
XIN)=-V*DXfW + DX 
YIM=-V«DY<-Z + CY 
YIM)=-V^OYfZ + DY 
RETURN 
ENC 

SUBROUTINE PLOTT I , CV ) 

COMMON/INTFAC/XI 1800 ) , YI 1800) 

LGGICAL-n MINUS,LABL 
COMMGN/TABL/ TA3C I 20 ,6 ) , JC 

COMMCjN/DITS/XMI.\',YMIN, SLOPEX, SLCPEYfOITSDX, DITSDY , 
SCALE POINTS FOR PLOT ROUTINE 
DO IOC 1=1 , NP 
XI I ) = SLOPEX^(XI I )-XMIN) 

YI I ) =SLGPE Y*(-YI I )-YMI N) 

CALL LINEIX ,Y,NP,1 , 1 ) 

IFI .NCT.LAbL) RETURN 
DIR=O.C 

GO TO 11,2,3,4,6), lOIR 
DIR=9C. 

CCCRDX=XII ) 

C0GRDY=Y(1 ) 

CALL NUMBER I COOR DX , C OORDY , . 07 , C V ,0 I R , 3 ) 

RETURN 

MOVE PEN DOWN ONE HALF INCH 
DIR=90. 

CGCR0X=XI1 ) 

CCCRDY=YIi )*.3 
GO TO 5 

MOVE PEN TO THE LEFT 
CCCRDX=XI1 )-.3 
C0CR0Y=Y(1 ) 

GO TO 5 

SEARCH FOR XMAX , XM I N , YMAX , YM I N , AND SAVE YMINX 
XMAX=XI 1 ) 

SMIN=XMAX 
YMIMX=YI1 ) 

YMAX=YMIKX 

VM IN = YMIi\X 

DO 200 1=2, NP 

IFIXI I) .GT.XMAX) XMAX = Xin 

IFIYI D.LT.VMIN) VMIN = vm 

IFIYI I ) .GT.YMAX) YMAX = YII) 

IFIXI I ) .GE. SHIN) GO TO 200 
SMIN=XII) 

YMINX=Y( I ) 

CONTINUE 

JC=NUMBER OF ENTRIES IN TABC 
IFIJC) 40C, 500,400 
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400 DO 900 1=1, JC 

IF(XMAX.LT.TABC( I ,U . AND . YMAX . L T ,TA BC < I , 2 ) . A ND. VM IN . 
900 CONTINUE 

DID NOT FIND THIS CONTOUR TO BE INTERIOR TO ANOTHER 
CHECK IF EXTERIOR 
DO iCCC I =1 , JC 

IFIXMAX.GT.TABCI 1,1) . AND . YMAX . GT . T ABC ( 1 , 2 J . A ND . VM IN . 
1000 CONTINUE 
500 IF (JC.EQ.20) RETURN 
JC=JC+1 
MC = JC 

600 TABC { MC, 1 ) =XMAX 
TAEC( MC,2 ) =YRAX 
TA0C( RC,3)=VNIN 
TABC(MC,4)=SMIN 
TABC(NC,5)=YMINX 
TA6C ( NC,6) =CV 
RETURN 

C CHECK IF THIS INTERIOR ONE IS OF HIGHER LEVEL 

700 IF(CV.LE.TA6C(I,6) ) RETURN 
2000 MC=I 

GO TO 600 

C CHECK IF LEVEL GF THIS EXTERIOR ONE IS HIGHER 

800 IF(CV.LT.TABC( 1,6) ) RETURN 
GO TO 2000 
END 
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3 * l2iiSiSEl Program 

The prograh which solves equation A. 18 from Appendix A 

is shown in the following pages. It calculates the solution, 

to the boundary layer voltage drop for a varying C . It 

P 

therefore requires an enthalpy table for the calculation, 
and such a table is provided in data format following the 
program. This particular table is for a 

Toluene/oxygen/cesiura mixture at temperatures from 1500<>K to 
26400K. An explanation on the use of the prograra is 
included . 
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c 

c 

c 



c 



BDUNOARY LAYFR ROUTINE 

THIS ROUTINE IS DESIGNED TO CALCULATE THE BOUNDARY 
LAYER VOLTAGE DROP IN AN MHD CHANNEL. IT USES A 
SIMPSON'S PULE INTEGRATION TECHNIQUE AND REQUIRES AN 
ENTHALPY TABLE FOR THE RANGE OF TEMPERATURES CONSIDERED, 
THE FCLLCWING PARAMETERS ARE REQUIRED: 

POT=SFED IONIZATION POTENTIAL IN EV'S 
TFRFE^FPEE STREAM TEMPERATURE 
TWALL=hALL TEMPERATURE 

TTCP=MAX TEMPERATURE IN ENTHALPY TABLE. TTOP MUST BE 
GREATER THAN OR EQUAL TO TFREE. 

TBCT=LCWEST TEMPERATURE IN ENTHALPY TABLE. MUST BE 
EQUAL TO OR SMALLER THAN TWALL. 

NPTS=i\UMBER OF POINTS USED IN INTEGRATION ROUTINE 
MAX=NURBER CF POINTS IN ENTHALPY TABLE 
ENTHALPY TABLE MUST BE EQUALLY INCREMENTED 

COMMOK/BLl/OEL , T FR E E , TNA LL , TBOT , AL PH , ENT Z , T T OP , MAX,N,E 
NTH( 100 ) 

READ (E , lOU POT, TFP E E , T-a'A LL ,NPT S , N 
WRITE(6, 100) POT, TFREE, TWALL, NPTStN 
READ(5,i02> TTCP,TBOT,MAX 
READ! 5,103) (ENTHI J) , J=--l ,MAX) 

WRITE (6, 105) 

MAXI =MAX-1 



NN = N-1 



NM=N-2 

ALPH=PCT-5. 8C54E 03/TFREE 
DEL= (TTOP-TBOT) /FLOAT (MAX- 1 ) 

THWALL=TWALL/TFPEE 
DELM=( TWALL- TBOT ) /DEL 
M=DELK 

ENTZ = F.NTH(f1 + l ) + ( ENTH ( M + 2 ) -ENTH ( M + 1 ) ) ^' ( TWALL- M*D EL - 
TBCT) /DEL 

DELM= { TFREE-TBOT )/DEL 
M=CELP 

IFIM.LT.MAXI ) GO TO 4 
ENTF = ENTH(MAX ) 

GO TO 5 

4 ENTF = ENTH(M + l )^■{ ENTH ( M + 2 ) - ENTH ( M+ 1 ) ) ’’^ ( TFR EE - M’^D EL- 



TBCT) /DEL 

5 ENTHM1=ENTF-ENTZ 
FAC = ENTHM1--*N 
CO EF 2 = FLO AT (N) -TFREE /FAC 
C=FCTM 1.0) 

DELTA=( 1 ,0-THWALL)/FL0aT( NPTS-1 ) 

SIMPSON'S RULE INTEGRATION PATTERN 1,4, 2, 4,1 
SUM1=0.0 
SUM2=C.O 
NPT1=NPTS-1 



NFT2=NFTS-2 

COMPUTATION FOR LAST POINT. 
SUM=C 

COMPUTATION FOR EVEN POINTS 
THET A=THWALL+DELTA 
DC 10 J=2,NPTi,2 
SUM 1 = SUM H-FCTN( THETA ) 
THET A=THETA +2 .-DELTA 
10 CONTINUE 

COMPUTATION FOR ODD POINTS. 
THETA = THWALL+-2.-DELTA 
DO 20 J=3,NPT2,2 
SUM2=SUM2+FCTN<THETA ) 
THETA=THETA+2.*DcLTA 



FIRST POINT IS ZERO. 
..THOSE MULTIPLIED BY 



.THOSE MULTIPLIED BY 



4 



2 



20 CONTINUE ^ , ,, 

VQLT = DELT A*( SUM+4.*SUMl + 2 .^^SUM2 )/3. 

VGLT=CCEr2*V0LT-l .0 

WR IT E (6, 104 J TWALL ,T FREE , THWALL , VOLT 
STOP 

LOl FCFMAT (3F9. C, 219 ) 

LOO FORMAT (3X?' IONIZATION POTENT I AL = * , F 7 ? 3 , 3X , ' F RF E STRF AM 
TEMP = ' ,F6. 1 , 3X, ' WALL TEMP = ' , F6 . 1 , iX , ' N P TS = ' , 14 , ' 
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103 FCPMAT (7E9. 0) 

104 F0RMAT{3X,2F9.0tF9.3,E12.4) 

10 5 FORM AT ( /,7X , • TWALL ' , 4X , • TFRFE' , 3X , • THWALL ' , 3X , ' VOLT* } 
END 
C 

C FUNCTION THETA 
C 

FUNCTION FCTN( THETA) 

COMMCN/BL 1 / DEL »T FREE , TWALL, TBOT,AL PH, ENTZtTTOP ,MAX,M,E 
NTHdOO) 

MAX1=MAX-1 

MAXM=KAX-3 

NN=M-1 

A = THETAS--* (-1.75) 

B=EXP( ALPH^M 1 .O/THFTA-1.0) ) 

TEFP = THETA=5=TFRFE 
DELM=^( TEMP-T30T) /DEL 
M=DELF 

I F( M.LT.MAXI ) go to 4 
ENT=ENTH(MAX ) 

GO TO 5 

4 ENT = ENTii( f-Hi ) + ( ENTH( M+^2 )-ENTH( M + 1 ) ) =!- ( TEMP-M«OE L - 

TBCT)/DEL 

5 ENT1=ENT-ENTZ 
C = ENTH-*NN 
IF(M.LT.2) GO TO 1 

I FCM.GT.MAX'v’) GO TO 2 

CP = ( ENTH( M- 1 )-6.^ENTH( M)+8.^:‘ENTH(M+-2 )-ENTH(M + 
3))/(12.^D£L) 

GO TO 3 

1 C? = { -50.*ENTH( M<-1 ) f 96. ^ENTH( M + 2 )-7 2.X=ENTH( iU3 ) + 

22.*FNTH(M+4)-6. *ENTH( M+ 

5))/(24.^^-DEL) 

GO TO 3 

2 CP = ( 6»«ENTH(M-5)-3 2. *ENTH( M-2) <-72.=!=EN7H(M-l ) - 

96.^ENTH( M) +50. =:=EfMTH( M+- 
1 ) )/( 24.’:=DEL} 

3 FCTN = A-S*C=:=CP 
RETURN 

END 



C DATA 

3.87 260C. 2000. 401 4 

2640. 15CC. 58 



-1745.28 


- 1738. 22 
17C2.58 


-1731.14 


-1724.04 


-1716.91 


-1709.76 - 


-1695. 36 


- 1688. 10 
1650.61 


-1680. 79 


-1673.41 


-1665.93 


-1658.34 - 


-1642. 7C 


-1634.59 

1590.33 


-1626.26 


-1617.69 


-1603.85 


-1599.73 - 


-1580.62 


-1570.59 
1515. 05 


-1560.23 


-1549.51 


-1538.42 


-1526.95 - 


-1502.72 


- 1489.93 
1418.11 


-1476.65 


-1462.85 


-1448.52 


-1433.61 - 


- 1401.98 


-1285, 18 
) 290.29 


-1367. 71 


-1349.51 


-1330.55 


-1310.83 - 


-1268.90 


-1246.64 
1 121. 20 


-1223.48 


-1199.33 


-1174. 32 


-1148.27 - 


-1093.09 

-864.93 


-1C63.90 

9C1.04 

-827.60 


-1033.61 


-1002.21 


-969.66 


-935.94 
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